Merge pull request #163 from tansongchen/main

Add Householder's method

Author: ChrisRackauckas

Project.toml

@@ -26,12 +26,14 @@ ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
 Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
+TaylorDiff = "b36ab563-344f-407b-a36a-4f200bebf99c"
 
 [extensions]
 SimpleNonlinearSolveChainRulesCoreExt = "ChainRulesCore"
 SimpleNonlinearSolveReverseDiffExt = "ReverseDiff"
 SimpleNonlinearSolveTrackerExt = "Tracker"
 SimpleNonlinearSolveZygoteExt = "Zygote"
+SimpleNonlinearSolveTaylorDiffExt = "TaylorDiff"
 
 [compat]
 ADTypes = "1.9"
@@ -63,6 +65,7 @@ SciMLBase = "2.37.0"
 Setfield = "1.1.1"
 StaticArrays = "1.9"
 StaticArraysCore = "1.4.2"
+TaylorDiff = "0.2.5"
 Test = "1.10"
 Tracker = "0.2.33"
 Zygote = "0.6.69"
@@ -86,9 +89,10 @@ ReTestItems = "817f1d60-ba6b-4fd5-9520-3cf149f6a823"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+TaylorDiff = "b36ab563-344f-407b-a36a-4f200bebf99c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["AllocCheck", "Aqua", "DiffEqBase", "ExplicitImports", "FiniteDiff", "ForwardDiff", "Hwloc", "InteractiveUtils", "LinearAlgebra",  "NonlinearProblemLibrary", "Pkg", "PolyesterForwardDiff", "Random", "ReTestItems", "Reexport", "ReverseDiff",  "StaticArrays", "Test", "Tracker", "Zygote"]
+test = ["AllocCheck", "Aqua", "DiffEqBase", "ExplicitImports", "FiniteDiff", "ForwardDiff", "Hwloc", "InteractiveUtils", "LinearAlgebra",  "NonlinearProblemLibrary", "Pkg", "PolyesterForwardDiff", "Random", "ReTestItems", "Reexport", "ReverseDiff",  "StaticArrays", "TaylorDiff", "Test", "Tracker", "Zygote"]

ext/SimpleNonlinearSolveTaylorDiffExt.jl

@@ -0,0 +1,62 @@
+module SimpleNonlinearSolveTaylorDiffExt
+using SimpleNonlinearSolve
+using SimpleNonlinearSolve: ImmutableNonlinearProblem, ReturnCode, build_solution,
+                            check_termination, init_termination_cache
+using SimpleNonlinearSolve: __maybe_unaliased, _get_fx, __fixed_parameter_function
+using MaybeInplace: @bb
+using SciMLBase: isinplace
+
+import TaylorDiff
+
+@inline function __get_higher_order_derivatives(
+        ::SimpleHouseholder{N}, prob, f, x, fx) where {N}
+    vN = Val(N)
+    l = map(one, x)
+    t = TaylorDiff.make_seed(x, l, vN)
+
+    if isinplace(prob)
+        bundle = similar(fx, TaylorDiff.TaylorScalar{eltype(fx), N})
+        f(bundle, t)
+        map!(TaylorDiff.primal, fx, bundle)
+    else
+        bundle = f(t)
+        fx = map(TaylorDiff.primal, bundle)
+    end
+    bundle = inv.(bundle)
+    num = TaylorDiff.extract_derivative(bundle, N - 1)
+    den = TaylorDiff.extract_derivative(bundle, N)
+    return num, den, fx
+end
+
+function SciMLBase.__solve(prob::ImmutableNonlinearProblem, alg::SimpleHouseholder{N},
+        args...; abstol = nothing, reltol = nothing, maxiters = 1000,
+        termination_condition = nothing, alias_u0 = false, kwargs...) where {N}
+    x = __maybe_unaliased(prob.u0, alias_u0)
+    length(x) == 1 ||
+        throw(ArgumentError("SimpleHouseholder only supports scalar problems"))
+    fx = _get_fx(prob, x)
+    @bb xo = copy(x)
+    f = __fixed_parameter_function(prob)
+
+    abstol, reltol, tc_cache = init_termination_cache(
+        prob, abstol, reltol, fx, x, termination_condition)
+
+    for i in 1:maxiters
+        num, den, fx = __get_higher_order_derivatives(alg, prob, f, x, fx)
+
+        if i == 1
+            iszero(fx) && build_solution(prob, alg, x, fx; retcode = ReturnCode.Success)
+        else
+            # Termination Checks
+            tc_sol = check_termination(tc_cache, fx, x, xo, prob, alg)
+            tc_sol !== nothing && return tc_sol
+        end
+
+        @bb copyto!(xo, x)
+        @bb x .+= (N - 1) .* num ./ den
+    end
+
+    return build_solution(prob, alg, x, fx; retcode = ReturnCode.MaxIters)
+end
+
+end

src/SimpleNonlinearSolve.jl

@@ -47,6 +47,7 @@ include("nlsolve/lbroyden.jl")
 include("nlsolve/klement.jl")
 include("nlsolve/trustRegion.jl")
 include("nlsolve/halley.jl")
+include("nlsolve/householder.jl")
 include("nlsolve/dfsane.jl")
 
 ## Interval Nonlinear Solvers
@@ -139,6 +140,7 @@ end
 export AutoFiniteDiff, AutoForwardDiff, AutoPolyesterForwardDiff
 export SimpleBroyden, SimpleDFSane, SimpleGaussNewton, SimpleHalley, SimpleKlement,
        SimpleLimitedMemoryBroyden, SimpleNewtonRaphson, SimpleTrustRegion
+export SimpleHouseholder
 export Alefeld, Bisection, Brent, Falsi, ITP, Ridder
 
 end # module

src/nlsolve/householder.jl

@@ -0,0 +1,16 @@
+"""
+    SimpleHouseholder{order}()
+
+A low-overhead implementation of Householder's method to arbitrary order.
+This method is non-allocating on scalar and static array problems.
+
+!!! warning
+
+    Needs `TaylorDiff.jl` to be explicitly loaded before using this functionality.
+    Internally, this uses TaylorDiff.jl for automatic differentiation.
+
+### Type Parameters
+
+  - `order`: the convergence order of the Householder method. `order = 2` is the same as Newton's method, `order = 3` is the same as Halley's method, etc.
+"""
+struct SimpleHouseholder{order} <: AbstractNewtonAlgorithm end

test/core/rootfind_tests.jl

@@ -1,6 +1,7 @@
 @testsetup module RootfindingTesting
 using Reexport
-@reexport using AllocCheck, StaticArrays, Random, LinearAlgebra, ForwardDiff, DiffEqBase
+@reexport using AllocCheck, StaticArrays, Random, LinearAlgebra, ForwardDiff, DiffEqBase,
+                TaylorDiff
 import PolyesterForwardDiff
 
 quadratic_f(u, p) = u .* u .- p
@@ -91,6 +92,32 @@ end
     end
 end
 
+@testitem "SimpleHouseholder" setup=[RootfindingTesting] tags=[:core] begin
+    @testset "AutoDiff: TaylorDiff.jl" for order in (2, 3, 4)
+        @testset "[OOP] u0: $(nameof(typeof(u0)))" for u0 in ([1.0], @SVector[1.0], 1.0)
+            sol = benchmark_nlsolve_oop(
+                quadratic_f, u0; solver = SimpleHouseholder{order}())
+            @test SciMLBase.successful_retcode(sol)
+            @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
+        end
+
+        @testset "[IIP] u0: $(nameof(typeof(u0)))" for u0 in ([1.0],)
+            sol = benchmark_nlsolve_iip(
+                quadratic_f!, u0; solver = SimpleHouseholder{order}())
+            @test SciMLBase.successful_retcode(sol)
+            @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
+        end
+    end
+
+    @testset "Termination condition: $(nameof(typeof(termination_condition))) u0: $(nameof(typeof(u0)))" for termination_condition in TERMINATION_CONDITIONS,
+        u0 in (1.0, [1.0], @SVector[1.0])
+
+        probN = NonlinearProblem(quadratic_f, u0, 2.0)
+        @test all(solve(probN, SimpleHouseholder{2}(); termination_condition).u .≈
+                  sqrt(2.0))
+    end
+end
+
 @testitem "Derivative Free Metods" setup=[RootfindingTesting] tags=[:core] begin
     @testset "$(nameof(typeof(alg)))" for alg in [
         SimpleBroyden(), SimpleKlement(), SimpleDFSane(),