Skip to content

WIP: Integration of Anderson acceleration of fixed-point iteration #12

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
wants to merge 1 commit into from
Closed

WIP: Integration of Anderson acceleration of fixed-point iteration #12

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 4 additions & 1 deletion src/DelayDiffEq.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ using OrdinaryDiffEq, DataStructures, RecursiveArrayTools, Combinatorics

using Compat

using Base.Test
using NLsolve

import OrdinaryDiffEq: initialize!, perform_step!, loopfooter!,
loopheader!, alg_order, handle_tstop!, ODEIntegrator, savevalues!,
Expand All @@ -22,6 +22,9 @@ import DiffEqBase: solve, solve!, init, resize!, u_cache, user_cache,
import OrdinaryDiffEq: Rosenbrock23Cache, Rosenbrock32Cache,
ImplicitEulerCache, TrapezoidCache

import NLsolve: AbstractDifferentiableMultivariateFunction

include("iteration_function.jl")
include("integrator_type.jl")
include("integrator_interface.jl")
include("history_function.jl")
Expand Down
20 changes: 11 additions & 9 deletions src/algorithms.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,15 +3,17 @@

immutable MethodOfSteps{algType,AType,RType,NType,constrained} <: AbstractMethodOfStepsAlgorithm{constrained}
alg::algType
picardabstol::AType
picardreltol::RType
picardnorm::NType
max_picard_iters::Int
fixedpoint_abstol::AType
fixedpoint_reltol::RType
max_fixedpoint_iters::Int
picardnorm::NType # anderson acceleration with nlsolve always uses infinite norm of residuals
m::Int # controls history size of anderson acceleration (m=0 corresponds to simple fixed-point iteration)
end

Base.@pure MethodOfSteps(alg;constrained=false,
picardabstol = nothing,
picardreltol = nothing,
picardnorm = nothing,
max_picard_iters = 10) =
MethodOfSteps{typeof(alg),typeof(picardabstol),typeof(picardreltol),typeof(picardnorm),constrained}(alg,picardabstol,picardreltol,picardnorm,max_picard_iters)
fixedpoint_abstol = nothing,
fixedpoint_reltol = nothing,
max_fixedpoint_iters = 10,
picardnorm = nothing,
m = 0) =
MethodOfSteps{typeof(alg),typeof(fixedpoint_abstol),typeof(fixedpoint_reltol),typeof(picardnorm),constrained}(alg,fixedpoint_abstol,fixedpoint_reltol,max_fixedpoint_iters,picardnorm,m)
99 changes: 67 additions & 32 deletions src/integrator_interface.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,67 +1,102 @@
function savevalues!(integrator::DDEIntegrator,force_save=false)
# update ODE integrator with values of DDE integrator
integrator.integrator.u = integrator.u
integrator.integrator.k = integrator.k
integrator.integrator.t = integrator.t

# add steps for interpolation when needed
OrdinaryDiffEq.ode_addsteps!(integrator.integrator,integrator.f)

# update solution of ODE integrator
savevalues!(integrator.integrator,force_save)
end

function postamble!(integrator::DDEIntegrator)
# update ODE integrator with values of DDE integrator
integrator.integrator.u = integrator.u
integrator.integrator.k = integrator.k
integrator.integrator.t = integrator.t

# clean up solution of ODE integrator
OrdinaryDiffEq.postamble!(integrator.integrator)
end

function perform_step!(integrator::DDEIntegrator)
integrator.tprev = integrator.t # this is necessary to extrapolate from the current interval

# update ODE integrator with values of DDE integrator
integrator.integrator.uprev = integrator.uprev
integrator.integrator.tprev = integrator.tprev
integrator.integrator.fsalfirst = integrator.fsalfirst
integrator.integrator.t = integrator.t
integrator.integrator.dt = integrator.dt

# if dt>max lag, then it's explicit so use Picard iteration
# if dt>min lag, then it's explicit so use fixed-point iteration
if integrator.dt >minimum(integrator.prob.lags)

# the is done to correct the extrapolation
t_prev_cache = integrator.tprev
# this is done to correct the extrapolation
t_cache = integrator.t
uprev_cache = integrator.uprev
tprev_cache = integrator.tprev
if typeof(integrator.uprev_cache) <: AbstractArray
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I thought it might be useful to save allocations for uprev_cache and add a new field to the integrator.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That doesn't save allocations since uprev_cache = integrator.uprev is just saving the pointer to another variable.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh yes, my bad. So the question is rather: why is it saved in the first place? If entries of integrator.uprev are changed, entries of uprev_cache are also changed - so when we later set integrator.uprev to uprev_cache again, entries of integrator.uprev are in fact not updated to any prior values. Moreover, just by skimming through low_order_rk_integrators.jl in OrdinaryDiffEq.jl it seems that integrator.uprev is not updated by perform_step! - but I assume there are algorithms which also update uprev? Otherwise uprev_cache could be removed completely.

copy!(integrator.uprev_cache,integrator.uprev)
else
integrator.uprev_cache = integrator.uprev
end

numiters = 1
while true
if typeof(integrator.u) <: AbstractArray
copy!(integrator.u_cache,integrator.u)
else
integrator.u_cache = integrator.u
end
perform_step!(integrator,integrator.cache)
if eltype(integrator.u) <: Real # not possible to use nlsolve with units (yet)
integrator.first_iteration = true # first iteration step needs special updates

if typeof(integrator.resid) <: AbstractArray
integrator.resid .= (integrator.u .- integrator.u_cache)./(integrator.picardabstol .+ max.(abs.(integrator.u),abs.(integrator.u_cache)).*integrator.picardreltol)
# execute Anderson acceleration of fixed-point iteration
if typeof(integrator.u) <: Vector
nlsolve(integrator.iterator,integrator.u;
method=:anderson,m=integrator.m,iterations=integrator.max_fixedpoint_iters,xtol=0,ftol=1)
else
integrator.resid = (integrator.u .- integrator.u_cache)./(integrator.picardabstol .+ max.(abs.(integrator.u),abs.(integrator.u_cache)).*integrator.picardreltol)
end

picardEEst = integrator.picardnorm(integrator.resid)
if picardEEst < 1 || numiters > integrator.max_picard_iters
break
nlsolve(integrator.iterator,vec(collect(integrator.u));
method=:anderson,m=integrator.m,iterations=integrator.max_fixedpoint_iters,xtol=0,ftol=1)
end
if numiters == 1
integrator.integrator.tprev = integrator.t
integrator.integrator.t = integrator.t+integrator.dt
integrator.integrator.uprev = integrator.u
else # simple Picard iteration for units
for i in 1:integrator.max_fixedpoint_iters
# update cache of u
if typeof(integrator.u) <: AbstractArray
copy!(integrator.u_cache,integrator.u)
else
integrator.u_cache = integrator.u
end

perform_step!(integrator,integrator.cache)

# calculate residuals
if typeof(integrator.resid) <: AbstractArray
@. integrator.resid = (integrator.u - integrator.u_cache)/
@muladd(integrator.fixedpoint_abstol + max(abs(integrator.u),abs(integrator.u_cache))*integrator.fixedpoint_reltol)
else
integrator.resid = @. (integrator.u - integrator.u_cache)/
(integrator.fixedpoint_abstol + max(abs(integrator.u),abs(integrator.u_cache))*integrator.fixedpoint_reltol)
end

# special updates of ODE integrator in first iteration step
if i == 1
integrator.integrator.tprev = integrator.t
integrator.integrator.t = integrator.t+integrator.dt
integrator.integrator.uprev = integrator.u
end

# update ODE integrator with values of DDE integrator
integrator.integrator.u = integrator.u
integrator.integrator.k = integrator.k

# stop fixed-point iteration when residuals are small
integrator.picardnorm(integrator.resid) < 1 && break
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Currently the exit condition of the loop is checked directly after calculating the residuals. However, then no updates would be applied to the ODE integrator if the loop is left during the first iteration? For now, I moved it to the end - the same way Anderson accelerated fixed-point iteration with NLsolve does it right now.

end
integrator.integrator.u = integrator.u
integrator.integrator.k = integrator.k
numiters += 1
end

# reset to cached values
integrator.t = t_cache
integrator.integrator.t = t_cache
integrator.tprev = t_prev_cache
integrator.uprev = uprev_cache

integrator.tprev = tprev_cache
if typeof(integrator.uprev) <: AbstractArray
copy!(integrator.uprev,integrator.uprev_cache)
else
integrator.uprev = integrator.uprev_cache
end
else # no iterations
perform_step!(integrator,integrator.cache)
end
Expand Down
27 changes: 18 additions & 9 deletions src/integrator_type.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,11 +9,13 @@ type DDEIntegrator{algType<:OrdinaryDiffEqAlgorithm,uType,tType,absType,relType,
uprev::uType
tprev::tType
u_cache::uType
picardabstol::absType
picardreltol::relType
uprev_cache::uType
fixedpoint_abstol::absType
fixedpoint_reltol::relType
resid::residType # This would have to resize for resizing DDE to work
picardnorm::NType
max_picard_iters::Int
max_fixedpoint_iters::Int
m::Int
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe a more meaningful should be used instead of m?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I prefer descriptive variable names. history_size?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That sounds good 👍

alg::algType
rate_prototype::rateType
notsaveat_idxs::Vector{Int}
Expand All @@ -39,26 +41,33 @@ type DDEIntegrator{algType<:OrdinaryDiffEqAlgorithm,uType,tType,absType,relType,
integrator::IType
fsalfirst::rateType
fsallast::rateType
first_iteration::Bool
iterator::IterationFunction{DDEIntegrator{algType,uType,tType,absType,relType,residType,tTypeNoUnits,tdirType,ksEltype,SolType,rateType,F,ProgressType,CacheType,IType,ProbType,NType,O}}

(::Type{DDEIntegrator{algType,uType,tType,absType,relType,
residType,tTypeNoUnits,tdirType,ksEltype,SolType,rateType,F,
ProgressType,CacheType,IType,ProbType,NType,O}}){
algType,uType,tType,absType,
relType,residType,tTypeNoUnits,tdirType,ksEltype,SolType,rateType,F,ProgressType,
CacheType,IType,ProbType,NType,O}(sol,prob,u,k,t,dt,f,uprev,tprev,u_cache,
picardabstol,picardreltol,resid,picardnorm,max_picard_iters,
CacheType,IType,ProbType,NType,O}(sol,prob,u,k,t,dt,f,uprev,tprev,u_cache,uprev_cache,
fixedpoint_abstol,fixedpoint_reltol,resid,picardnorm,max_fixedpoint_iters,m,
alg,rate_prototype,notsaveat_idxs,dtcache,dtchangeable,dtpropose,tdir,
EEst,qold,q11,
iter,saveiter,saveiter_dense,prog,cache,
kshortsize,just_hit_tstop,accept_step,isout,reeval_fsal,u_modified,
integrator,opts) = new{algType,uType,tType,absType,relType,residType,tTypeNoUnits,tdirType,
opts,integrator) = begin
dde_int = new{algType,uType,tType,absType,relType,residType,tTypeNoUnits,tdirType,
ksEltype,SolType,rateType,F,ProgressType,CacheType,IType,ProbType,NType,O}(
sol,prob,u,k,t,dt,f,uprev,tprev,u_cache,
picardabstol,picardreltol,resid,picardnorm,max_picard_iters,
sol,prob,u,k,t,dt,f,uprev,tprev,u_cache,uprev_cache,
fixedpoint_abstol,fixedpoint_reltol,resid,picardnorm,max_fixedpoint_iters,m,
alg,rate_prototype,notsaveat_idxs,dtcache,dtchangeable,dtpropose,tdir,
EEst,qold,q11,
iter,saveiter,saveiter_dense,prog,cache,
kshortsize,just_hit_tstop,accept_step,isout,reeval_fsal,u_modified,integrator,opts) # Leave off fsalfirst and last
kshortsize,just_hit_tstop,accept_step,isout,reeval_fsal,u_modified,opts,integrator) # Leave off fsalfirst, fsallast, first_iteration and iterator
dde_int.iterator = IterationFunction{DDEIntegrator{algType,uType,tType,absType,relType,residType,tTypeNoUnits,tdirType,
ksEltype,SolType,rateType,F,ProgressType,CacheType,IType,ProbType,NType,O}}(dde_int)
dde_int
end
end

function (integrator::DDEIntegrator)(t,deriv::Type=Val{0};idxs=nothing)
Expand Down
41 changes: 41 additions & 0 deletions src/iteration_function.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,41 @@
# nlsolve optimizes subtypes of AbstractDifferentiableMultivariateFunction
mutable struct IterationFunction{I<:AbstractDDEIntegrator} <: AbstractDifferentiableMultivariateFunction
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

"Regular" functions are always encapsulated DifferentiableMultivariateFunction by NLsolve. I'm not sure whether there is an advantage in defining a subtype of AbstractDifferentiableMultivariateFunction with the self-referential field f! instead of just subtyping Function

integrator::I
f!::IterationFunction{I} # Anderson acceleration needs field f!

function IterationFunction{I}(integrator::I) where I <:AbstractDDEIntegrator
f = new(integrator)
f.f! = f
end
end

function (f::IterationFunction)(x,fvec)
integrator = f.integrator

# update u
if typeof(integrator.u) <: AbstractArray
copy!(integrator.u,x)
else
integrator.u = x[1] # x is always a vector
end

perform_step!(integrator,integrator.cache)

# update output vector of residuals
# in contrast to picard iteration we can not use integrator.resid
# as output vector for different fixed-point iterations
@. fvec = (integrator.u - x)/
@muladd(integrator.fixedpoint_abstol + max(abs(integrator.u),abs(x))*integrator.fixedpoint_reltol)

# special updates of ODE integrator after the first iteration necessary
if integrator.first_iteration
integrator.integrator.tprev = integrator.t
integrator.integrator.t = integrator.t+integrator.dt
integrator.integrator.uprev = integrator.u
integrator.first_iteration = false
end

# update ODE integrator with values of DDE integrator
integrator.integrator.u = integrator.u
integrator.integrator.k = integrator.k
end
37 changes: 23 additions & 14 deletions src/solve.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -73,10 +73,10 @@ function init{uType,tType,isinplace,algType<:AbstractMethodOfStepsAlgorithm,lTyp
end
integrator.dt = dt

if typeof(alg.picardabstol) <: Void
picardabstol_internal = map(eltype(uType),integrator.opts.abstol)
if typeof(alg.fixedpoint_abstol) <: Void
fixedpoint_abstol_internal = map(eltype(uType),integrator.opts.abstol)
else
picardabstol_internal = map(eltype(uType),alg.picardabstol)
fixedpoint_abstol_internal = map(eltype(uType),alg.fixedpoint_abstol)
end
if typeof(alg.picardnorm) <: Void
picardnorm = integrator.opts.internalnorm
Expand All @@ -85,23 +85,32 @@ function init{uType,tType,isinplace,algType<:AbstractMethodOfStepsAlgorithm,lTyp

uEltypeNoUnits = typeof(recursive_one(integrator.u))

if typeof(alg.picardreltol) <: Void
picardreltol_internal = map(uEltypeNoUnits,integrator.opts.reltol)
if typeof(alg.fixedpoint_reltol) <: Void
fixedpoint_reltol_internal = map(uEltypeNoUnits,integrator.opts.reltol)
else
picardreltol_internal = map(uEltypeNoUnits,alg.picardreltol)
fixedpoint_reltol_internal = map(uEltypeNoUnits,alg.fixedpoint_reltol)
end
if typeof(integrator.u) <: AbstractArray
resid = similar(integrator.u,uEltypeNoUnits)
u_cache = similar(integrator.u)
uprev_cache = similar(integrator.u)
else
resid = one(uEltypeNoUnits)
u_cache = oneunit(eltype(uType))
uprev_cache = oneunit(eltype(uType))
end

# for real numbers nlsolve is used for Anderson acceleration of fixed-point iteration which creates vectors of residuals
if eltype(integrator.u) <: Real
resid = nothing
elseif typeof(integrator.u) <: AbstractArray
resid = similar(integrator.u,uEltypeNoUnits)
else
resid = one(uEltypeNoUnits)
end

dde_int = DDEIntegrator{typeof(integrator.alg),
uType,tType,
typeof(picardabstol_internal),
typeof(picardreltol_internal),
typeof(fixedpoint_abstol_internal),
typeof(fixedpoint_reltol_internal),
typeof(resid),
tTypeNoUnits,typeof(integrator.tdir),
typeof(integrator.k),typeof(sol),
Expand All @@ -112,17 +121,17 @@ function init{uType,tType,isinplace,algType<:AbstractMethodOfStepsAlgorithm,lTyp
typeof(picardnorm),
typeof(integrator.opts)}(
sol,prob,integrator.u,integrator.k,integrator.t,integrator.dt,
dde_f2,integrator.uprev,integrator.tprev,u_cache,
picardabstol_internal,picardreltol_internal,
resid,picardnorm,alg.max_picard_iters,
dde_f2,integrator.uprev,integrator.tprev,u_cache,uprev_cache,
fixedpoint_abstol_internal,fixedpoint_reltol_internal,
resid,picardnorm,alg.max_fixedpoint_iters,alg.m,
integrator.alg,integrator.rate_prototype,integrator.notsaveat_idxs,integrator.dtcache,
integrator.dtchangeable,integrator.dtpropose,integrator.tdir,
integrator.EEst,integrator.qold,integrator.q11,
integrator.iter,integrator.saveiter,
integrator.saveiter_dense,integrator.prog,integrator.cache,
integrator.kshortsize,integrator.just_hit_tstop,integrator.accept_step,
integrator.isout,
integrator.reeval_fsal,integrator.u_modified,integrator.opts,integrator) # Leave off fsalfirst and last
integrator.reeval_fsal,integrator.u_modified,integrator.opts,integrator) # Leave off fsalfirst, fasllast, first_iteration, and iterator

initialize!(dde_int)
initialize!(integrator.opts.callback,integrator.t,integrator.u,dde_int)
Expand Down
Loading