Add Householder's method

Author: tansongchen

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

lib/SimpleNonlinearSolve/ext/SimpleNonlinearSolveTaylorDiffExt.jl

@@ -0,0 +1,60 @@
+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) where N
+    vN = Val(N)
+    l = map(one, x)
+
+    if isinplace(prob)
+        fx = f(x)
+        invf = x -> inv.(f(x))
+        bundle = TaylorDiff.derivatives(invf, fx, x, l, vN)
+    else
+        t = TaylorDiff.make_seed(x, l, vN)
+        ft = f(t)
+        fx = map(TaylorDiff.primal, ft)
+        bundle = inv.(ft)
+    end
+    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)
+
+        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

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

lib/SimpleNonlinearSolve/src/nlsolve/householder.jl

@@ -0,0 +1,13 @@
+"""
+    SimpleHouseholder{order}()
+
+A low-overhead implementation of Householder's method. This method is non-allocating on scalar
+and static array problems.
+
+Internally, this uses TaylorDiff.jl for the 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