WIP: Integration of Anderson acceleration of fixed-point iteration #12
Conversation
| integrator.integrator.k = integrator.k | ||
|
|
||
| # stop fixed-point iteration when residuals are small | ||
| integrator.picardnorm(integrator.resid) < 1 && break |
There was a problem hiding this comment.
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.
| t_cache = integrator.t | ||
| uprev_cache = integrator.uprev | ||
| tprev_cache = integrator.tprev | ||
| if typeof(integrator.uprev_cache) <: AbstractArray |
There was a problem hiding this comment.
I thought it might be useful to save allocations for uprev_cache and add a new field to the integrator.
There was a problem hiding this comment.
That doesn't save allocations since uprev_cache = integrator.uprev is just saving the pointer to another variable.
There was a problem hiding this comment.
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.
| picardnorm::NType | ||
| max_picard_iters::Int | ||
| max_fixedpoint_iters::Int | ||
| m::Int |
There was a problem hiding this comment.
Maybe a more meaningful should be used instead of m?
There was a problem hiding this comment.
Yes, I prefer descriptive variable names. history_size?
|
Would it be too much to separate out the naming changes + comments to a separate PR? Those are really good, and I would like those merged immediately. And separating those will make reviewing maintaining the change to Andersson much easier (it'll be much easier to find out if anything in the future breaks due to that change if it's isolated). Of course, naming and comments which are only for Anderson should be in the PR with Anderson, but I see a lot of auxiliary changes which are riding along and making it a bit harder to follow. |
| @@ -0,0 +1,41 @@ | |||
| # nlsolve optimizes subtypes of AbstractDifferentiableMultivariateFunction | |||
| mutable struct IterationFunction{I<:AbstractDDEIntegrator} <: AbstractDifferentiableMultivariateFunction | |||
There was a problem hiding this comment.
"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
Codecov Report
@@ Coverage Diff @@
## master #12 +/- ##
==========================================
+ Coverage 79.26% 81.34% +2.07%
==========================================
Files 8 9 +1
Lines 164 193 +29
==========================================
+ Hits 130 157 +27
- Misses 34 36 +2
Continue to review full report at Codecov.
|
1 similar comment
|
I'll move the comments + naming changes to a seperate PR. Since I'm on a trip until the beginning of next week it might take some days, but hopefully I can do it tomorrow. |
|
Take your time. Have a nice trip. I think what's left here is resolving the conflicts and then the review. |
|
I would like to complete #17 first. If it is merged it is probably easier to update this PR, and also the changes can be restricted to the implementation of Anderson acceleration only. |
|
I wonder if this is the right approach. NLsolve seems to allocate a bunch: https://github.com/JuliaNLSolvers/NLsolve.jl/blob/master/src/anderson.jl#L17 and it'll do a dynamic dispatch each time too. This makes it fine for longer fixed point iteration problems, but for quick problems that's going to have a big overhead. Given that an the fact that most of the changes in this PR are due to commenting which was another PR, I think this can be closed. However, one approach to getting Anderson acceleration may be to rebase this PR, get it done, then test a new Anderson acceleration implementation against it. |
|
It would be preferable to allocate these arrays once in the beginning of our algorithm, so they could be reused in every fixed-point iteration (although there is still a not in-place |
|
But I also guess we should close this PR for now. |
Yeah, I feel bad for asking for it in NLsolve but yes it seems we need a few specialized changes so that will make their implementation not work for our usage. |
4 participants
This is a first draft trying to integrate Anderson acceleration. Unfortunately (currently) NLsolve accepts only vectors of real numbers, so I had to keep the simple Picard iteration for compatibility with Unitful.jl. Locally the tests passed already, however it is not benchmarked yet. I also tried to improve readability by adding some comments.