Started to implement DFSane

Author: Deltadahl

src/SimpleNonlinearSolve.jl

@@ -24,13 +24,14 @@ include("klement.jl")
 include("trustRegion.jl")
 include("ridder.jl")
 include("brent.jl")
+include("dfsane.jl")
 include("ad.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, Broyden, Klement, SimpleTrustRegion)
+    for alg in (SimpleNewtonRaphson, Broyden, Klement, SimpleTrustRegion, SimpleDFSane)
         solve(prob_no_brack, alg(), abstol = T(1e-2))
     end
 
@@ -51,7 +52,7 @@ SnoopPrecompile.@precompile_all_calls begin for T in (Float32, Float64)
 end end
 
 # DiffEq styled algorithms
-export Bisection, Brent, Broyden, Falsi, Klement, Ridder, SimpleNewtonRaphson,
+export Bisection, Brent, Broyden, SimpleDFSane, Falsi, Klement, Ridder, SimpleNewtonRaphson,
        SimpleTrustRegion
 
 end # module

src/ad.jl

@@ -31,16 +31,16 @@ end
 function SciMLBase.solve(prob::NonlinearProblem{<:Union{Number, StaticArraysCore.SVector},
                                                 iip,
                                                 <:Dual{T, V, P}},
-                         alg::Union{SimpleNewtonRaphson, SimpleTrustRegion},
+                         alg::Union{SimpleNewtonRaphson, SimpleTrustRegion, SimpleDFSane}, # TODO make this to one variable
                          args...; kwargs...) where {iip, T, V, P}
     sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...)
     return SciMLBase.build_solution(prob, alg, Dual{T, V, P}(sol.u, partials), sol.resid;
                                     retcode = sol.retcode)
 end
-function SciMLBase.solve(prob::NonlinearProblem{<:Union{Number, StaticArraysCore.SVector},
+function SciMLBase.solve(prob::NonlinearProblem{<:Union{Number, StaticArraysCore.SVector}, # TODO make this to one variable
                                                 iip,
                                                 <:AbstractArray{<:Dual{T, V, P}}},
-                         alg::Union{SimpleNewtonRaphson, SimpleTrustRegion}, args...;
+                         alg::Union{SimpleNewtonRaphson, SimpleTrustRegion, SimpleDFSane}, args...;
                          kwargs...) where {iip, T, V, P}
     sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...)
     return SciMLBase.build_solution(prob, alg, Dual{T, V, P}(sol.u, partials), sol.resid;

src/dfsane.jl

@@ -0,0 +1,135 @@
+"""
+```julia
+SimpleDFSane(; σ_min::Real = 1e-10, σ_0::Real = 1.0, M::Int = 10,
+            γ::Real = 1e-4, τ_min::Real = 0.1, τ_max::Real = 0.5,
+            nexp::Int = 2)
+```
+
+A low-overhead implementation of the df-sane method. For more information, see [1].
+References:
+    W LaCruz, JM Martinez, and M Raydan (2006), Spectral residual mathod without gradient information for solving large-scale nonlinear systems of equations, Mathematics of Computation, 75, 1429-1448.
+### Keyword Arguments
+
+
+- `σ_min`: the minimum value of `σ_k`. Defaults to `1e-10`. # TODO write about this...
+- `σ_0`: the initial value of `σ_k`. Defaults to `1.0`. # TODO write about this...
+- `M`: the value of `M` in the paper. Defaults to `10`. # TODO write about this...
+- `γ`: the value of `γ` in the paper. Defaults to `1e-4`. # TODO write about this...
+- `τ_min`: the minimum value of `τ_k`. Defaults to `0.1`. # TODO write about this...
+- `τ_max`: the maximum value of `τ_k`. Defaults to `0.5`. # TODO write about this...
+- `nexp`: the value of `nexp` in the paper. Defaults to `2`. # TODO write about this...
+
+"""
+struct SimpleDFSane{T} <: AbstractSimpleNonlinearSolveAlgorithm
+    σ_min::T
+    σ_0::T
+    M::Int
+    γ::T
+    τ_min::T
+    τ_max::T
+    nexp::Int
+
+    function SimpleDFSane(; σ_min::Real = 1e-10, σ_0::Real = 1.0, M::Int = 10,
+                          γ::Real = 1e-4, τ_min::Real = 0.1, τ_max::Real = 0.5,
+                          nexp::Int = 2)
+        new{typeof(σ_min)}(σ_min, σ_0, M, γ, τ_min, τ_max, nexp)
+    end
+end
+
+function SciMLBase.__solve(prob::NonlinearProblem,
+                           alg::SimpleDFSane, args...; abstol = nothing,
+                           reltol = nothing,
+                           maxiters = 1000, kwargs...)
+    f = Base.Fix2(prob.f, prob.p)
+    x = float(prob.u0)
+    T = eltype(x)
+    σ_min = float(alg.σ_min)
+    σ_k = float(alg.σ_0)
+    M = alg.M
+    γ = float(alg.γ)
+    τ_min = float(alg.τ_min)
+    τ_max = float(alg.τ_max)
+    nexp = alg.nexp
+
+    if SciMLBase.isinplace(prob)
+        error("SimpleDFSane 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)
+
+    function ff(x)
+        F = f(x)
+        f_k = norm(F)^nexp
+        return f_k, F
+    end
+
+    f_k, F_k = ff(x)
+    α_0 = convert(T, 1.0)
+    f_0 = f_k
+    prev_fs = fill(f_k, M)
+
+    for k in 1:maxiters
+        iszero(F_k) &&
+            return SciMLBase.build_solution(prob, alg, x, F_k;
+                                            retcode = ReturnCode.Success)
+
+        # Control spectral parameter
+        if abs(σ_k) > 1 / σ_min
+            σ_k = 1 / σ_min * sign(σ_k)
+        elseif abs(σ_k) < σ_min
+            σ_k = σ_min
+        end
+
+        # Line search direction
+        d = -σ_k * F_k
+
+        # Nonmonotone line search
+        η = f_0 / k^2
+
+        f_bar = maximum(prev_fs)
+        α_p = α_0
+        α_m = α_0
+        xp = x + α_p * d
+        fp, Fp = ff(xp)
+        while true
+            if fp ≤ f_bar + η - γ * α_p^2 * f_k
+                break
+            end
+
+            α_tp = α_p^2 * f_k / (fp + (2 * α_p - 1) * f_k)
+            xp = x - α_m * d
+            fp, Fp = ff(xp)
+
+            if fp ≤ f_bar + η - γ * α_m^2 * f_k
+                break
+            end
+
+            α_tm = α_m^2 * f_k / (fp + (2 * α_m - 1) * f_k)
+            α_p = min(τ_max * α_p, max(α_tp, τ_min * α_p))
+            α_m = min(τ_max * α_m, max(α_tm, τ_min * α_m))
+            xp = x + α_p * d
+            fp, Fp = ff(xp)
+        end
+
+        if isapprox(xp, x, atol = atol, rtol = rtol)
+            return SciMLBase.build_solution(prob, alg, xp, Fp;
+                                            retcode = ReturnCode.Success)
+        end
+        # Update spectral parameter
+        s_k = xp - x
+        y_k = Fp - F_k
+        σ_k = dot(s_k, s_k) / dot(s_k, y_k)
+
+        # Take step
+        x = xp
+        F_k = Fp
+        f_k = fp
+
+        # Store function value
+        idx_to_replace = k % M + 1
+        prev_fs[idx_to_replace] = fp
+    end
+    return SciMLBase.build_solution(prob, alg, x, F_k; retcode = ReturnCode.MaxIters)
+end

test/basictests.jl

@@ -62,37 +62,49 @@ sol = benchmark_scalar(sf, csu0)
 @test sol.retcode === ReturnCode.Success
 @test sol.u * sol.u - 2 < 1e-9
 
+# SimpleDFSane
+function benchmark_scalar(f, u0)
+    probN = NonlinearProblem{false}(f, u0)
+    sol = (solve(probN, SimpleDFSane()))
+end
+
+sol = benchmark_scalar(sf, csu0)
+@test sol.retcode === ReturnCode.Success
+@test sol.u * sol.u - 2 < 1e-9
+
 # AD Tests
 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(), Klement(), SimpleTrustRegion(),
+            SimpleDFSane())
     g = function (p)
         probN = NonlinearProblem{false}(f, csu0, p)
         sol = solve(probN, alg, abstol = 1e-9)
         return sol.u[end]
     end
 
     for p in 1.1:0.1:100.0
-        @test g(p) ≈ sqrt(p)
-        @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p))
+        @test abs.(g(p)) ≈ sqrt(p)
+        @test abs.(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(), Klement(), SimpleTrustRegion())
+for alg in (SimpleNewtonRaphson(), Broyden(), Klement(), SimpleTrustRegion(),
+            SimpleDFSane())
     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)
-        @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p))
+        @test abs(g(p)) ≈ sqrt(p)
+        @test abs(ForwardDiff.derivative(g, p)) ≈ 1 / (2 * sqrt(p))
     end
 end
 
@@ -148,7 +160,7 @@ 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(), Klement(), SimpleTrustRegion(), SimpleDFSane())
     global g, p
     g = function (p)
         probN = NonlinearProblem{false}(f, 0.5, p)
@@ -169,6 +181,7 @@ probN = NonlinearProblem(f, u0)
 @test solve(probN, SimpleTrustRegion(; autodiff = false)).u[end] ≈ sqrt(2.0)
 @test solve(probN, Broyden()).u[end] ≈ sqrt(2.0)
 @test solve(probN, Klement()).u[end] ≈ sqrt(2.0)
+@test solve(probN, SimpleDFSane()).u[end] ≈ sqrt(2.0)
 
 for u0 in [1.0, [1, 1.0]]
     local f, probN, sol
@@ -185,8 +198,8 @@ 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, Klement()).u ≈ sol
+    @test solve(probN, SimpleDFSane()).u ≈ sol
 end
 
 # Bisection Tests
@@ -299,3 +312,46 @@ for options in list_of_options
     sol = solve(probN, alg)
     @test all(abs.(f(u, p)) .< 1e-10)
 end
+
+# # Test that `SimpleDFSane` passes a test that `SimpleNewtonRaphson` fails on.
+# u0 = [-10.0, -1.0, 1.0, 2.0, 3.0, 4.0, 10.0]
+# global g, f
+# f = (u, p) -> 0.010000000000000002 .+
+#               10.000000000000002 ./ (1 .+
+#                (0.21640425613334457 .+
+#                 216.40425613334457 ./ (1 .+
+#                  (0.21640425613334457 .+
+#                   216.40425613334457 ./
+#                   (1 .+ 0.0006250000000000001(u .^ 2.0))) .^ 2.0)) .^ 2.0) .-
+#               0.0011552453009332421u .- p
+# g = function (p)
+#     probN = NonlinearProblem{false}(f, u0, p)
+#     sol = solve(probN, SimpleDFSane())
+#     return sol.u
+# end
+# p = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+# u = g(p)
+# f(u, p)
+# @test all(abs.(f(u, p)) .< 1e-10)
+
+# # Test kwars in `SimpleDFSane`
+
+
+# list_of_options = zip(max_trust_radius, initial_trust_radius, step_threshold,
+#                       shrink_threshold, expand_threshold, shrink_factor,
+#                       expand_factor, max_shrink_times)
+# for options in list_of_options
+#     local probN, sol, alg
+#     alg = SimpleDFSane(max_trust_radius = options[1],
+#                             initial_trust_radius = options[2],
+#                             step_threshold = options[3],
+#                             shrink_threshold = options[4],
+#                             expand_threshold = options[5],
+#                             shrink_factor = options[6],
+#                             expand_factor = options[7],
+#                             max_shrink_times = options[8])
+
+#     probN = NonlinearProblem(f, u0, p)
+#     sol = solve(probN, alg)
+#     @test all(abs.(f(u, p)) .< 1e-10)
+# end