Skip to content

Commit

Permalink
Merge pull request #337 from avik-pal/ap/fixed_point
Browse files Browse the repository at this point in the history
FixedPointSolvers: Recommendations and Wrappers
  • Loading branch information
ChrisRackauckas authored Dec 23, 2023
2 parents d723578 + 29917d8 commit 1ada921
Show file tree
Hide file tree
Showing 15 changed files with 440 additions and 17 deletions.
12 changes: 10 additions & 2 deletions Project.toml
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "NonlinearSolve"
uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
authors = ["SciML"]
version = "3.1.2"
version = "3.2.0"

[deps]
ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
Expand Down Expand Up @@ -32,18 +32,22 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
[weakdeps]
BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
FastLevenbergMarquardt = "7a0df574-e128-4d35-8cbd-3d84502bf7ce"
FixedPointAcceleration = "817d07cb-a79a-5c30-9a31-890123675176"
LeastSquaresOptim = "0fc2ff8b-aaa3-5acd-a817-1944a5e08891"
MINPACK = "4854310b-de5a-5eb6-a2a5-c1dee2bd17f9"
NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
SpeedMapping = "f1835b91-879b-4a3f-a438-e4baacf14412"
Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"

[extensions]
NonlinearSolveBandedMatricesExt = "BandedMatrices"
NonlinearSolveFastLevenbergMarquardtExt = "FastLevenbergMarquardt"
NonlinearSolveFixedPointAccelerationExt = "FixedPointAcceleration"
NonlinearSolveLeastSquaresOptimExt = "LeastSquaresOptim"
NonlinearSolveMINPACKExt = "MINPACK"
NonlinearSolveNLsolveExt = "NLsolve"
NonlinearSolveSpeedMappingExt = "SpeedMapping"
NonlinearSolveSymbolicsExt = "Symbolics"
NonlinearSolveZygoteExt = "Zygote"

Expand All @@ -60,6 +64,7 @@ Enzyme = "0.11.11"
FastBroadcast = "0.2.8"
FastLevenbergMarquardt = "0.1"
FiniteDiff = "2.21"
FixedPointAcceleration = "0.3"
ForwardDiff = "0.10.36"
LazyArrays = "1.8.2"
LeastSquaresOptim = "0.8.5"
Expand All @@ -84,6 +89,7 @@ SciMLOperators = "0.3.7"
SimpleNonlinearSolve = "1.0.2"
SparseArrays = "<0.0.1, 1"
SparseDiffTools = "2.14"
SpeedMapping = "0.3"
StableRNGs = "1"
StaticArrays = "1.7"
Symbolics = "5.13"
Expand All @@ -99,6 +105,7 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
FastLevenbergMarquardt = "7a0df574-e128-4d35-8cbd-3d84502bf7ce"
FixedPointAcceleration = "817d07cb-a79a-5c30-9a31-890123675176"
ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
LeastSquaresOptim = "0fc2ff8b-aaa3-5acd-a817-1944a5e08891"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Expand All @@ -112,11 +119,12 @@ Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
SpeedMapping = "f1835b91-879b-4a3f-a438-e4baacf14412"
StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"

[targets]
test = ["Aqua", "Enzyme", "BenchmarkTools", "SafeTestsets", "Pkg", "Test", "ForwardDiff", "StaticArrays", "Symbolics", "LinearSolve", "Random", "LinearAlgebra", "Zygote", "SparseDiffTools", "NonlinearProblemLibrary", "LeastSquaresOptim", "FastLevenbergMarquardt", "NaNMath", "BandedMatrices", "DiffEqBase", "StableRNGs", "MINPACK", "NLsolve", "OrdinaryDiffEq"]
test = ["Aqua", "Enzyme", "BenchmarkTools", "SafeTestsets", "Pkg", "Test", "ForwardDiff", "StaticArrays", "Symbolics", "LinearSolve", "Random", "LinearAlgebra", "Zygote", "SparseDiffTools", "NonlinearProblemLibrary", "LeastSquaresOptim", "FastLevenbergMarquardt", "NaNMath", "BandedMatrices", "DiffEqBase", "StableRNGs", "MINPACK", "NLsolve", "OrdinaryDiffEq", "SpeedMapping", "FixedPointAcceleration"]
1 change: 1 addition & 0 deletions docs/pages.jl
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@ pages = ["index.md",
"solvers/BracketingSolvers.md",
"solvers/SteadyStateSolvers.md",
"solvers/NonlinearLeastSquaresSolvers.md",
"solvers/FixedPointSolvers.md",
"solvers/LineSearch.md"],
"Detailed Solver APIs" => Any["api/nonlinearsolve.md",
"api/simplenonlinearsolve.md",
Expand Down
17 changes: 17 additions & 0 deletions docs/src/api/fixedpointacceleration.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,17 @@
# FixedPointAcceleration.jl

This is a extension for importing solvers from FixedPointAcceleration.jl into the SciML
interface. Note that these solvers do not come by default, and thus one needs to install
the package before using these solvers:

```julia
using Pkg
Pkg.add("FixedPointAcceleration")
using FixedPointAcceleration, NonlinearSolve
```

## Solver API

```@docs
FixedPointAccelerationJL
```
17 changes: 17 additions & 0 deletions docs/src/api/speedmapping.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,17 @@
# SpeedMapping.jl

This is a extension for importing solvers from SpeedMapping.jl into the SciML
interface. Note that these solvers do not come by default, and thus one needs to install
the package before using these solvers:

```julia
using Pkg
Pkg.add("SpeedMapping")
using SpeedMapping, NonlinearSolve
```

## Solver API

```@docs
SpeedMappingJL
```
42 changes: 42 additions & 0 deletions docs/src/solvers/FixedPointSolvers.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,42 @@
# Fixed Point Solvers

Currently we don't have an API to directly specify Fixed Point Solvers. However, a Fixed
Point Problem can be trivially converted to a Root Finding Problem. Say we want to solve:

```math
f(u) = u
```

This can be written as:

```math
g(u) = f(u) - u = 0
```

``g(u) = 0`` is a root finding problem. Note that we can use any root finding
algorithm to solve this problem. However, this is often not the most efficient way to
solve a fixed point problem. We provide a few algorithms available via extensions that
are more efficient for fixed point problems.

Note that even if you use one of the Fixed Point Solvers mentioned here, you must still
use the `NonlinearProblem` API to specify the problem, i.e., ``g(u) = 0``.

## Recommended Methods

Using [native NonlinearSolve.jl methods](@ref nonlinearsystemsolvers) is the recommended
approach. For systems where constructing Jacobian Matrices are expensive, we recommend
using a Krylov Method with one of those solvers.

## Full List of Methods

We are only listing the methods that natively solve fixed point problems.

### SpeedMapping.jl

- `SpeedMappingJL()`: accelerates the convergence of a mapping to a fixed point by the
Alternating cyclic extrapolation algorithm (ACX).

### FixedPointAcceleration.jl

- `FixedPointAccelerationJL()`: accelerates the convergence of a mapping to a fixed point
by the Anderson acceleration algorithm and a few other methods.
2 changes: 1 addition & 1 deletion docs/src/solvers/NonlinearSystemSolvers.md
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
# [Nonlinear System Solvers](@id nonlinearsystemsolvers)

`solve(prob::NonlinearProblem,alg;kwargs)`
`solve(prob::NonlinearProblem, alg; kwargs)`

Solves for ``f(u)=0`` in the problem defined by `prob` using the algorithm
`alg`. If no algorithm is given, a default algorithm will be chosen.
Expand Down
56 changes: 56 additions & 0 deletions ext/NonlinearSolveFixedPointAccelerationExt.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,56 @@
module NonlinearSolveFixedPointAccelerationExt

using NonlinearSolve, FixedPointAcceleration, DiffEqBase, SciMLBase

function SciMLBase.__solve(prob::NonlinearProblem, alg::FixedPointAccelerationJL, args...;
abstol = nothing, maxiters = 1000, alias_u0::Bool = false,
show_trace::Val{PrintReports} = Val(false), termination_condition = nothing,
kwargs...) where {PrintReports}
@assert (termination_condition ===
nothing)||(termination_condition isa AbsNormTerminationMode) "SpeedMappingJL does not support termination conditions!"

u0 = NonlinearSolve.__maybe_unaliased(prob.u0, alias_u0)
u_size = size(u0)
T = eltype(u0)
iip = isinplace(prob)
p = prob.p

if !iip && prob.u0 isa Number
# FixedPointAcceleration makes the scalar problem into a vector problem
f = (u) -> [prob.f(u[1], p) .+ u[1]]
elseif !iip && prob.u0 isa AbstractVector
f = (u) -> (prob.f(u, p) .+ u)
elseif !iip && prob.u0 isa AbstractArray
f = (u) -> vec(prob.f(reshape(u, u_size), p) .+ u)
elseif iip && prob.u0 isa AbstractVector
du = similar(u0)
f = (u) -> (prob.f(du, u, p); du .+ u)
else
du = similar(u0)
f = (u) -> (prob.f(du, reshape(u, u_size), p); vec(du) .+ u)
end

tol = abstol === nothing ? real(oneunit(T)) * (eps(real(one(T))))^(4 // 5) : abstol

sol = fixed_point(f, NonlinearSolve._vec(u0); Algorithm = alg.algorithm,
ConvergenceMetricThreshold = tol, MaxIter = maxiters, MaxM = alg.m,
ExtrapolationPeriod = alg.extrapolation_period, Dampening = alg.dampening,
PrintReports, ReplaceInvalids = alg.replace_invalids,
ConditionNumberThreshold = alg.condition_number_threshold, quiet_errors = true)

res = prob.u0 isa Number ? first(sol.FixedPoint_) : sol.FixedPoint_
if res === missing
resid = NonlinearSolve.evaluate_f(prob, u0)
res = u0
converged = false
else
resid = NonlinearSolve.evaluate_f(prob, res)
converged = maximum(abs, resid) tol
end
return SciMLBase.build_solution(prob, alg, res, resid;
retcode = converged ? ReturnCode.Success : ReturnCode.Failure,
stats = SciMLBase.NLStats(sol.Iterations_, 0, 0, 0, sol.Iterations_),
original = sol)
end

end
13 changes: 8 additions & 5 deletions ext/NonlinearSolveMINPACKExt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@ using MINPACK

function SciMLBase.__solve(prob::Union{NonlinearProblem{uType, iip},
NonlinearLeastSquaresProblem{uType, iip}}, alg::CMINPACK, args...;
abstol = 1e-6, maxiters = 100000, alias_u0::Bool = false,
abstol = nothing, maxiters = 100000, alias_u0::Bool = false,
termination_condition = nothing, kwargs...) where {uType, iip}
@assert (termination_condition ===
nothing)||(termination_condition isa AbsNormTerminationMode) "CMINPACK does not support termination conditions!"
Expand All @@ -16,6 +16,7 @@ function SciMLBase.__solve(prob::Union{NonlinearProblem{uType, iip},
u0 = NonlinearSolve.__maybe_unaliased(prob.u0, alias_u0)
end

T = eltype(u0)
sizeu = size(prob.u0)
p = prob.p

Expand All @@ -25,11 +26,11 @@ function SciMLBase.__solve(prob::Union{NonlinearProblem{uType, iip},

if !iip && prob.u0 isa Number
f! = (du, u) -> (du .= prob.f(first(u), p); Cint(0))
elseif !iip && prob.u0 isa Vector{Float64}
elseif !iip && prob.u0 isa AbstractVector
f! = (du, u) -> (du .= prob.f(u, p); Cint(0))
elseif !iip && prob.u0 isa AbstractArray
f! = (du, u) -> (du .= vec(prob.f(reshape(u, sizeu), p)); Cint(0))
elseif prob.u0 isa Vector{Float64}
elseif prob.u0 isa AbstractVector
f! = (du, u) -> prob.f(du, u, p)
else # Then it's an in-place function on an abstract array
f! = (du, u) -> (prob.f(reshape(du, sizeu), reshape(u, sizeu), p); du = vec(du); 0)
Expand All @@ -43,14 +44,16 @@ function SciMLBase.__solve(prob::Union{NonlinearProblem{uType, iip},
method = ifelse(alg.method === :auto,
ifelse(prob isa NonlinearLeastSquaresProblem, :lm, :hybr), alg.method)

abstol = abstol === nothing ? real(oneunit(T)) * (eps(real(one(T))))^(4 // 5) : abstol

if SciMLBase.has_jac(prob.f)
if !iip && prob.u0 isa Number
g! = (du, u) -> (du .= prob.f.jac(first(u), p); Cint(0))
elseif !iip && prob.u0 isa Vector{Float64}
elseif !iip && prob.u0 isa AbstractVector
g! = (du, u) -> (du .= prob.f.jac(u, p); Cint(0))
elseif !iip && prob.u0 isa AbstractArray
g! = (du, u) -> (du .= vec(prob.f.jac(reshape(u, sizeu), p)); Cint(0))
elseif prob.u0 isa Vector{Float64}
elseif prob.u0 isa AbstractVector
g! = (du, u) -> prob.f.jac(du, u, p)
else # Then it's an in-place function on an abstract array
g! = function (du, u)
Expand Down
16 changes: 10 additions & 6 deletions ext/NonlinearSolveNLsolveExt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,8 +3,9 @@ module NonlinearSolveNLsolveExt
using NonlinearSolve, NLsolve, DiffEqBase, SciMLBase
import UnPack: @unpack

function SciMLBase.__solve(prob::NonlinearProblem, alg::NLsolveJL, args...; abstol = 1e-6,
maxiters = 1000, alias_u0::Bool = false, termination_condition = nothing, kwargs...)
function SciMLBase.__solve(prob::NonlinearProblem, alg::NLsolveJL, args...;
abstol = nothing, maxiters = 1000, alias_u0::Bool = false,
termination_condition = nothing, kwargs...)
@assert (termination_condition ===
nothing)||(termination_condition isa AbsNormTerminationMode) "NLsolveJL does not support termination conditions!"

Expand All @@ -14,6 +15,7 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::NLsolveJL, args...; abst
u0 = NonlinearSolve.__maybe_unaliased(prob.u0, alias_u0)
end

T = eltype(u0)
iip = isinplace(prob)

sizeu = size(prob.u0)
Expand All @@ -25,11 +27,11 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::NLsolveJL, args...; abst

if !iip && prob.u0 isa Number
f! = (du, u) -> (du .= prob.f(first(u), p); Cint(0))
elseif !iip && prob.u0 isa Vector{Float64}
elseif !iip && prob.u0 isa AbstractVector
f! = (du, u) -> (du .= prob.f(u, p); Cint(0))
elseif !iip && prob.u0 isa AbstractArray
f! = (du, u) -> (du .= vec(prob.f(reshape(u, sizeu), p)); Cint(0))
elseif prob.u0 isa Vector{Float64}
elseif prob.u0 isa AbstractVector
f! = (du, u) -> prob.f(du, u, p)
else # Then it's an in-place function on an abstract array
f! = (du, u) -> (prob.f(reshape(du, sizeu), reshape(u, sizeu), p); du = vec(du); 0)
Expand All @@ -46,11 +48,11 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::NLsolveJL, args...; abst
if SciMLBase.has_jac(prob.f)
if !iip && prob.u0 isa Number
g! = (du, u) -> (du .= prob.f.jac(first(u), p); Cint(0))
elseif !iip && prob.u0 isa Vector{Float64}
elseif !iip && prob.u0 isa AbstractVector
g! = (du, u) -> (du .= prob.f.jac(u, p); Cint(0))
elseif !iip && prob.u0 isa AbstractArray
g! = (du, u) -> (du .= vec(prob.f.jac(reshape(u, sizeu), p)); Cint(0))
elseif prob.u0 isa Vector{Float64}
elseif prob.u0 isa AbstractVector
g! = (du, u) -> prob.f.jac(du, u, p)
else # Then it's an in-place function on an abstract array
g! = function (du, u)
Expand All @@ -68,6 +70,8 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::NLsolveJL, args...; abst
df = OnceDifferentiable(f!, vec(u0), vec(resid); autodiff)
end

abstol = abstol === nothing ? real(oneunit(T)) * (eps(real(one(T))))^(4 // 5) : abstol

original = nlsolve(df, vec(u0); ftol = abstol, iterations = maxiters, method,
store_trace, extended_trace, linesearch, linsolve, factor, autoscale, m, beta,
show_trace)
Expand Down
42 changes: 42 additions & 0 deletions ext/NonlinearSolveSpeedMappingExt.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,42 @@
module NonlinearSolveSpeedMappingExt

using NonlinearSolve, SpeedMapping, DiffEqBase, SciMLBase

function SciMLBase.__solve(prob::NonlinearProblem, alg::SpeedMappingJL, args...;
abstol = nothing, maxiters = 1000, alias_u0::Bool = false,
store_trace::Val{store_info} = Val(false), termination_condition = nothing,
kwargs...) where {store_info}
@assert (termination_condition ===
nothing)||(termination_condition isa AbsNormTerminationMode) "SpeedMappingJL does not support termination conditions!"

if typeof(prob.u0) <: Number
u0 = [prob.u0]
else
u0 = NonlinearSolve.__maybe_unaliased(prob.u0, alias_u0)
end

T = eltype(u0)
iip = isinplace(prob)
p = prob.p

if !iip && prob.u0 isa Number
m! = (du, u) -> (du .= prob.f(first(u), p) .+ first(u))
elseif !iip
m! = (du, u) -> (du .= prob.f(u, p) .+ u)
else
m! = (du, u) -> (prob.f(du, u, p); du .+= u)
end

tol = abstol === nothing ? real(oneunit(T)) * (eps(real(one(T))))^(4 // 5) : abstol

sol = speedmapping(u0; m!, tol, Lp = Inf, maps_limit = maxiters, alg.orders,
alg.check_obj, store_info, alg.σ_min, alg.stabilize)
res = prob.u0 isa Number ? first(sol.minimizer) : sol.minimizer
resid = NonlinearSolve.evaluate_f(prob, sol.minimizer)

return SciMLBase.build_solution(prob, alg, res, resid;
retcode = sol.converged ? ReturnCode.Success : ReturnCode.Failure,
stats = SciMLBase.NLStats(sol.maps, 0, 0, 0, sol.maps), original = sol)
end

end
3 changes: 2 additions & 1 deletion src/NonlinearSolve.jl
Original file line number Diff line number Diff line change
Expand Up @@ -236,7 +236,8 @@ export RadiusUpdateSchemes

export NewtonRaphson, TrustRegion, LevenbergMarquardt, DFSane, GaussNewton, PseudoTransient,
Broyden, Klement, LimitedMemoryBroyden
export LeastSquaresOptimJL, FastLevenbergMarquardtJL, CMINPACK, NLsolveJL
export LeastSquaresOptimJL,
FastLevenbergMarquardtJL, CMINPACK, NLsolveJL, FixedPointAccelerationJL, SpeedMappingJL
export NonlinearSolvePolyAlgorithm,
RobustMultiNewton, FastShortcutNonlinearPolyalg, FastShortcutNLLSPolyalg

Expand Down
Loading

0 comments on commit 1ada921

Please sign in to comment.