Merge pull request #344 from ErikQQY/qqy/siam_fixedpoint

Add SIAMFANLEquations fixed point solvers

Author: ChrisRackauckas

Project.toml

@@ -83,7 +83,7 @@ NonlinearProblemLibrary = "0.1.2"
 OrdinaryDiffEq = "6.63"
 Pkg = "1"
 PrecompileTools = "1.2"
-Printf = "1.9"
+Printf = "1.10"
 Random = "1.91"
 RecursiveArrayTools = "3.2"
 Reexport = "1.2"
@@ -92,7 +92,7 @@ SafeTestsets = "0.1"
 SciMLBase = "2.11"
 SciMLOperators = "0.3.7"
 SimpleNonlinearSolve = "1.0.2"
-SparseArrays = "1.9"
+SparseArrays = "1.10"
 SparseDiffTools = "2.14"
 SpeedMapping = "0.3"
 StableRNGs = "1"

docs/src/solvers/FixedPointSolvers.md

@@ -47,3 +47,7 @@ In our tests, we have found the anderson method implemented here to NOT be the m
 robust.
 
   - `NLsolveJL(; method = :anderson)`: Anderson acceleration for fixed point problems.
+
+### SIAMFANLEquations.jl
+
+  - `SIAMFANLEquationsJL(; method = :anderson)`: Anderson acceleration for fixed point problems.

ext/NonlinearSolveSIAMFANLEquationsExt.jl

@@ -24,7 +24,7 @@ end
 # pseudo transient continuation has a fixed cost per iteration, iteration statistics are
 # not interesting here.
 @inline function __siam_fanl_equations_stats_mapping(method, sol)
-    method === :pseudotransient && return nothing
+    ((method === :pseudotransient) || (method === :anderson)) && return nothing
     return SciMLBase.NLStats(sum(sol.stats.ifun), sum(sol.stats.ijac), 0, 0,
         sum(sol.stats.iarm))
 end
@@ -36,16 +36,14 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::SIAMFANLEquationsJL, arg
     @assert (termination_condition ===
              nothing)||(termination_condition isa AbsNormTerminationMode) "SIAMFANLEquationsJL does not support termination conditions!"
 
-    (; method, delta, linsolve) = alg
-
-    iip = SciMLBase.isinplace(prob)
+    (; method, delta, linsolve, m, beta) = alg
 
     T = eltype(prob.u0)
     atol = NonlinearSolve.DEFAULT_TOLERANCE(abstol, T)
     rtol = NonlinearSolve.DEFAULT_TOLERANCE(reltol, T)
 
     if prob.u0 isa Number
-        f = (u) -> prob.f(u, prob.p)
+        f = method == :anderson ? (du, u) -> (du = prob.f(u, prob.p)) : ((u) -> prob.f(u, prob.p))
 
         if method == :newton
             sol = nsolsc(f, prob.u0; maxit = maxiters, atol, rtol, printerr = ShT)
@@ -54,11 +52,16 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::SIAMFANLEquationsJL, arg
                 printerr = ShT)
         elseif method == :secant
             sol = secant(f, prob.u0; maxit = maxiters, atol, rtol, printerr = ShT)
+        elseif method == :anderson
+            f, u = NonlinearSolve.__construct_f(prob; alias_u0,
+                make_fixed_point = Val(true), can_handle_arbitrary_dims = Val(true))
+            sol = aasol(f, [prob.u0], m, __zeros_like(u, 1, 2*m+4); maxit = maxiters,
+                atol, rtol, beta = beta)
         end
 
         retcode = __siam_fanl_equations_retcode_mapping(sol)
         stats = __siam_fanl_equations_stats_mapping(method, sol)
-        resid = NonlinearSolve.evaluate_f(prob, sol.solution)
+        resid = NonlinearSolve.evaluate_f(prob, sol.solution[1])
         return SciMLBase.build_solution(prob, alg, sol.solution, resid; retcode,
             stats, original = sol)
     end
@@ -104,6 +107,11 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::SIAMFANLEquationsJL, arg
         elseif method == :pseudotransient
             sol = ptcsol(f!, u, FS, FPS; atol, rtol, maxit = maxiters,
                 delta0 = delta, printerr = ShT)
+        elseif method == :anderson
+            f!, u = NonlinearSolve.__construct_f(prob; alias_u0,
+                can_handle_arbitrary_dims = Val(true), make_fixed_point = Val(true))
+            sol = aasol(f!, u, m, zeros(T, N, 2*m+4), atol = atol, rtol = rtol,
+                maxit = maxiters, beta = beta)
         end
     else
         AJ!(J, u, x) = prob.f.jac(J, x, prob.p)

src/extension_algs.jl

@@ -377,23 +377,30 @@ end
   - `method`: the choice of method for solving the nonlinear system.
   - `delta`: initial pseudo time step, default is 1e-3.
   - `linsolve` : JFNK linear solvers, choices are `gmres` and `bicgstab`.
+  - `m`: Depth for Anderson acceleration, default as 0 for Picard iteration.
+  - `beta`: Anderson mixing parameter, change f(x) to (1-beta)x+beta*f(x),
+    equivalent to accelerating damped Picard iteration.
 
 ### Submethod Choice
 
   - `:newton`: Classical Newton method.
   - `:pseudotransient`: Pseudo transient method.
   - `:secant`: Secant method for scalar equations.
+  - `:anderson`: Anderson acceleration for fixed point iterations.
 """
 @concrete struct SIAMFANLEquationsJL{L <: Union{Symbol, Nothing}} <:
                  AbstractNonlinearSolveAlgorithm
     method::Symbol
     delta
     linsolve::L
+    m::Int
+    beta
 end
 
-function SIAMFANLEquationsJL(; method = :newton, delta = 1e-3, linsolve = nothing)
+function SIAMFANLEquationsJL(; method = :newton, delta = 1e-3, linsolve = nothing, m = 0,
+        beta = 1.0)
     if Base.get_extension(@__MODULE__, :NonlinearSolveSIAMFANLEquationsExt) === nothing
         error("SIAMFANLEquationsJL requires SIAMFANLEquations.jl to be loaded")
     end
-    return SIAMFANLEquationsJL(method, delta, linsolve)
+    return SIAMFANLEquationsJL(method, delta, linsolve, m, beta)
 end

test/wrappers/fixedpoint.jl

@@ -1,4 +1,4 @@
-using NonlinearSolve, FixedPointAcceleration, SpeedMapping, NLsolve, LinearAlgebra, Test
+using NonlinearSolve, FixedPointAcceleration, SpeedMapping, NLsolve, SIAMFANLEquations, LinearAlgebra, Test
 
 # Simple Scalar Problem
 @testset "Simple Scalar Problem" begin
@@ -14,6 +14,7 @@ using NonlinearSolve, FixedPointAcceleration, SpeedMapping, NLsolve, LinearAlgeb
     @test abs(solve(prob, SpeedMappingJL(; stabilize = true)).resid) ≤ 1e-10
 
     @test abs(solve(prob, NLsolveJL(; method = :anderson)).resid) ≤ 1e-10
+    @test abs(solve(prob, SIAMFANLEquationsJL(; method = :anderson)).resid) ≤ 1e-10
 end
 
 # Simple Vector Problem
@@ -28,6 +29,7 @@ end
     @test maximum(abs.(solve(prob, SpeedMappingJL()).resid)) ≤ 1e-10
     @test maximum(abs.(solve(prob, SpeedMappingJL(; orders = [3, 2])).resid)) ≤ 1e-10
     @test maximum(abs.(solve(prob, SpeedMappingJL(; stabilize = true)).resid)) ≤ 1e-10
+    @test maximum(abs.(solve(prob, SIAMFANLEquationsJL(; method = :anderson)).resid)) ≤ 1e-10
 
     @test_broken maximum(abs.(solve(prob, NLsolveJL(; method = :anderson)).resid)) ≤ 1e-10
 end
@@ -67,6 +69,9 @@ end
     sol = solve(prob, NLsolveJL(; method = :anderson))
     @test_broken sol.u' * A[:, 3] ≈ 32.916472867168096
 
+    sol = solve(prob, SIAMFANLEquationsJL(; method = :anderson))
+    @test sol.u' * A[:, 3] ≈ 32.916472867168096
+
     # Non vector inputs
     function power_method_nonvec!(du, u, A)
         mul!(vec(du), A, vec(u))