Predictor Corrector method #53

onoderat · 2018-02-14T05:04:01Z

The predictor corrector method for SDEs have been implemented and tested.

Both the in place and not in place versions have been implemented.
The testing code used the problem from the commutative_tests.jl .

The speed of the code has not been benchmarked rigorously. For a single setting of dt=0.002, the following speeds were achieved.

    EM in place            : 1.665ms (21771 allocations: 2.01 MiB)
    RKMil in place         : 1.881ms (23776 allocations: 2.10 MiB)
    PCEuler not in place   : 2.626ms (40673 allocations: 3.44 MiB)
    PCEuler in place       : 1.942ms (24271 allocations: 2.16 MiB)

The convergence tests shows that both the PC has better accuracy than EM and surprisingly also RKMilCommute. Although the PC is proven to be order 0.5, it showed order 1 scaling for this particular problem (see figure below). At dt=1/2^10, the l2 errors are respectively
RKMil : 6.02e-3
PCEuler : 1.44e-3

Added PCEuler as a type. alg_order and alg_compatible has been set. Has not set any other property - this may need to change.

The code has been refactored from basic_method_caches to the new file cache_types.jl. Code passed the default tests.

It is a singleton object since there are no static constants for this algorithm. Plan is to first focus on getting the ConstantCache version working, before creating the more optimized Cache version.

Now it's fields contains the bbprime functions. The optional arguments of theta and eta are set to 0.5 by default.

Code current runs for not in place problems. It has been checked to agree with EM. Careful testing of the deviation from the analytics is the next step - to check for higher order bugs.

Both is in place and not in place versions have been implemented. Test code is in test/predcorr_test.jl. This current version of script includes all the plotting and the timing. The next commit which will remove these plotting and timing features for brievity and runspeed. In terms of speed, here are some output from the test script. EM in place : 1.665ms (21771 allocations: 2.01 MiB) RKMil in place : 1.881ms (23776 allocations: 2.10 MiB) PC not in place: 2.626ms (40673 allocations: 3.44 MiB) PC in place : 1.942ms (24271 allocations: 2.16 MiB) The convergence test shows that both the PC has order 1 scaling in practice. In fact, the overall offset point is also lower for PC. At dt=1/2^10 RKMil : 6.02e-3 PC : 1.44e-3

The testing code has now been incoorporated into test/runtest.jl . Due to some conflicts in the use of f and g, the test script for predictor correcter now uses f_nondiag, g_nondiag. Pkg.test on StochasticDiffEq passed. A pull request will be submitted with this commit.

ChrisRackauckas · 2018-02-14T05:56:55Z

src/algorithms.jl

+struct PCEuler{T<:Real} <: StochasticDiffEqAlgorithm
+  theta::T
+  eta::T
+  bbprime::Function


Parameterize this for stability purposes.

Could you elaborate on this a little more? Would you like the PCEuler to be parametrized by the type of the function bbprime?

In Julia, every function is a different type and Function is just an abstract type. This allows specialization on functions and then allows things like inlining. If you ::Function in a type, then it just knows it's a <:Function and cannot optimize directly on it. So yup, just make that ::F and a free-roaming F parameter.

I see, will make the change.

ChrisRackauckas · 2018-02-14T05:57:11Z

src/algorithms.jl

+Reference: Stochastics and Dynamics, Vol. 8, No. 3 (2008) 561–581
+Note that the original paper has a typo in the definition of bbprime...
+"""
+# TODO: Check if this can be implemented with @pure.


This doesn't need @pure. It'll infer without it.

Thanks - I will remove the comment.

ChrisRackauckas · 2018-02-14T05:58:11Z

src/perform_step/predcorr.jl

+  end
+
+  f_n = f(uprev,p,t)
+  ubar = @muladd uprev + dt*f_n + gdW_n


@muladd on the function makes everything inside @muladd, so these aren't needed (there are probably some other functions in this library that need cleaning for that)

Awesome - I did note that adding this did not change the speed. Will remove them.

ChrisRackauckas · 2018-02-14T05:58:52Z

src/perform_step/predcorr.jl

+    gdW_np1 = g(ubar,p,tnp1).*dW
+  end
+
+  u .= @muladd uprev + (theta*dt)*fbar_np1 + ((1-theta)*dt)*fbar_n + eta*gdW_np1 + (1-eta)*gdW_n


You don't want to .= for numbers and static arrays. It'll work but be slower.

Sure, I will copy what you did for EM then.

ChrisRackauckas · 2018-02-14T06:01:04Z

test/predcorr_tests.jl

@@ -0,0 +1,89 @@
+using StochasticDiffEq, DiffEqProblemLibrary, DiffEqDevTools, Base.Test


Just rename this file as nondiagonal_tests.jl: I assume other algorithms will get added to it later.

ChrisRackauckas · 2018-02-14T06:03:32Z

Looks great! Just a few minor changes and I'll merge. I'm happily surprised it does that much better on this problem: that's a very big improvement! I wonder if it's order 1 on commutative noise problems or something like that: that could make sense.

ChrisRackauckas · 2018-02-14T06:14:28Z

src/algorithms.jl

+Note that the original paper has a typo in the definition of bbprime...
+"""
+# TODO: Check if this can be implemented with @pure.
+PCEuler(bbprime::Function; theta::T=1/2, eta::T=1/2) where {T<:Real} =


Don't ::Function in a type signature. It has a non-intuitive anti-specialization behavior, but also not all functions <:Function.

Wow, I have been doing that for a while now (in code I write), thanks for the tip! Will change.

In personal codes it's fine. In package codes:

struct MyF end (::MyF)(t) = 2t f = MyF()

f(t) works but f is not <:Function. You run into these weird cases the more usage things get. But if you don't do silly things yourself, you're safe 👍 (except for the non-specialization thing).

Ah that reminds me that I saw you did something like that to define the analytic solution to f. Ok will remove the ::Function in the type signature.

codecov · 2018-02-14T06:17:40Z

Codecov Report

❗ No coverage uploaded for pull request base (master@38fe1ca). Click here to learn what that means.
The diff coverage is 87.93%.

@@            Coverage Diff            @@
##             master      #53   +/-   ##
=========================================
  Coverage          ?   76.98%           
=========================================
  Files             ?       32           
  Lines             ?     3389           
  Branches          ?        0           
=========================================
  Hits              ?     2609           
  Misses            ?      780           
  Partials          ?        0

Impacted Files	Coverage Δ
src/caches/basic_method_caches.jl	`73.21% <ø> (ø)`
src/algorithms.jl	`96.96% <100%> (ø)`
src/alg_utils.jl	`94.38% <100%> (ø)`
src/caches/cache_types.jl	`100% <100%> (ø)`
src/perform_step/predcorr.jl	`85.36% <85.36%> (ø)`
src/caches/predcorr_caches.jl	`90% <90%> (ø)`

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 38fe1ca...6ba3629. Read the comment docs.

ChrisRackauckas · 2018-02-14T06:48:47Z

Just noticed bbprime should change to ggprime since the rest of the docs and functions refer to the diffusion term by g.

ChrisRackauckas · 2018-02-14T12:25:31Z

Good stuff, thanks. After seeing that from RKMil, I wonder how much of that is from the RK handling of the derivative handling. Would you mind doing a PR to do ggprime implementations of Milstein?

ChrisRackauckas · 2018-02-14T12:43:35Z

This needs to get added to the docs page http://docs.juliadiffeq.org/latest/solvers/sde_solve.html . It should probably be recommended when ggprime is known, and the definition of ggprime will need to go somewhere there for now. Eventually we'll want to handle ggprime like in:

http://docs.juliadiffeq.org/latest/features/performance_overloads.html

but that should wait on the change to DiffEqFunction where we can then have specific setups that are different for the different problem types.

onoderat · 2018-02-14T19:01:35Z

Hi Chris - I would like to do the ggprime implementation of Milstein. The issue is that the one in the textbook does not have implicit noise and neither do any papers explicitly construct this.

I am a physicist so I cant really do the proofs for these methods. Perhaps, you want to have a stab at the proof? Once I have the formula, I can do a PR for the method.

ChrisRackauckas · 2018-02-14T19:05:27Z

What do you mean by implicit noise? This shouldn't have any implicit terms.

onoderat · 2018-02-14T19:08:19Z

Ah sorry - I meant implicitness in the diffusion. Basically the one in the textbook does not have a eta parameter - it only has implicitness in the drift. See last sentence of screenshot below -

ChrisRackauckas · 2018-02-14T19:30:13Z

I was thinking of doing the simpler one first:

onoderat · 2018-02-14T19:48:26Z

Ah I see. Sure that should be easy to code up. But currently I am kind of rushing to prepare for the APS march meeting. Perhaps we could discuss this further when we meet in person next month?

I am thinking of working on the Wiktorsson approximations and creating a better interface for that stuff through NoiseProcess. Once that is done, perhaps we can multiple dispatch on the type of problem (commutative, diagonal, etc...).

This way we do not need to define RKMIlCommute, RKMil etc. What do you think?

ChrisRackauckas · 2018-02-14T19:53:51Z

I am thinking of working on the Wiktorsson approximations and creating a better interface for that stuff through NoiseProcess. Once that is done, perhaps we can multiple dispatch on the type of problem (commutative, diagonal, etc...). This way we do not need to define RKMIlCommute, RKMil etc. What do you think?

Yeah, share your idea when you can. As always, there's no rush. Could you open an issue to discuss this? I'm very curious about your idea here. I want to do something for this, and we should have a more constrained version of jumps that allow for these kinds of schemes from Platen.

ChrisRackauckas · 2018-02-14T19:54:48Z

BTW, schemes are usually drift-implicit because the inverse of the normal has infinite moments making it hard to prove convergence. There's things like truncated schemes though to get around this.

onoderat · 2018-02-14T20:08:16Z

Sure I will open an issue on the NoiseProcess stuff.

Thanks for explaining why the diffusion-implicit is hard. Interesting for the predictor-corrector, I numerically found that eta=0.5 is critical for getting that order 1 scaling.

Perhaps its an interesting research direction to investigate whether the same higher order scaling can be done if we have an eta (implicit diffusion) for the milstein case? Like say we get order 1.5 for doing order 1 level work?

ChrisRackauckas · 2018-02-14T20:10:43Z

That's super convergence. A lot of it is just finding the right conditions to make convergence "better", like how on additive noise the Euler-Maruyama scheme is Order 1.0. It would be very interesting to find out super convergence criterias for high order SDE schemes. I would assume that many of the common categories (additive, multiplicative) have some super convergence stuff going on.

onoderat added 7 commits February 9, 2018 12:59

Added predictor-corrector-euler algorithm type

42e2510

Added PCEuler as a type. alg_order and alg_compatible has been set. Has not set any other property - this may need to change.

Refactored definitions of general cache types

d27f23f

The code has been refactored from basic_method_caches to the new file cache_types.jl. Code passed the default tests.

Created ConstantCache type for predictor corrector alg.

71aa202

It is a singleton object since there are no static constants for this algorithm. Plan is to first focus on getting the ConstantCache version working, before creating the more optimized Cache version.

Modifield Predictor Corrector algorithm type

672f8a6

Now it's fields contains the bbprime functions. The optional arguments of theta and eta are set to 0.5 by default.

First version of predictor-corrector perform step

7872c2f

Code current runs for not in place problems. It has been checked to agree with EM. Careful testing of the deviation from the analytics is the next step - to check for higher order bugs.

ChrisRackauckas reviewed Feb 14, 2018

View reviewed changes

Made edits that were suggested in pull request.

6ba3629

ChrisRackauckas merged commit ccc5ca4 into SciML:master Feb 14, 2018

ChrisRackauckas mentioned this pull request Feb 14, 2018

Automatic ggprime handling #54

Open

ChrisRackauckas mentioned this pull request Feb 14, 2018

Implementing predictor corrector method for SDEs SciML/DifferentialEquations.jl#247

Closed

onoderat deleted the pull-request/ca05c51d branch February 14, 2018 20:16

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Predictor Corrector method #53

Predictor Corrector method #53

onoderat commented Feb 14, 2018 •

edited

ChrisRackauckas Feb 14, 2018

onoderat Feb 14, 2018 •

edited

ChrisRackauckas Feb 14, 2018

onoderat Feb 14, 2018

ChrisRackauckas Feb 14, 2018

onoderat Feb 14, 2018

ChrisRackauckas Feb 14, 2018

onoderat Feb 14, 2018

ChrisRackauckas Feb 14, 2018

onoderat Feb 14, 2018

ChrisRackauckas Feb 14, 2018

onoderat Feb 14, 2018

ChrisRackauckas commented Feb 14, 2018

ChrisRackauckas Feb 14, 2018

onoderat Feb 14, 2018 •

edited

ChrisRackauckas Feb 14, 2018

onoderat Feb 14, 2018

codecov bot commented Feb 14, 2018 •

edited

ChrisRackauckas commented Feb 14, 2018

ChrisRackauckas commented Feb 14, 2018

ChrisRackauckas commented Feb 14, 2018

onoderat commented Feb 14, 2018

ChrisRackauckas commented Feb 14, 2018

onoderat commented Feb 14, 2018

ChrisRackauckas commented Feb 14, 2018

onoderat commented Feb 14, 2018

ChrisRackauckas commented Feb 14, 2018

ChrisRackauckas commented Feb 14, 2018

onoderat commented Feb 14, 2018 •

edited

ChrisRackauckas commented Feb 14, 2018

		@@ -0,0 +1,89 @@
		using StochasticDiffEq, DiffEqProblemLibrary, DiffEqDevTools, Base.Test

Predictor Corrector method #53

Predictor Corrector method #53

Conversation

onoderat commented Feb 14, 2018 • edited

Choose a reason for hiding this comment

onoderat Feb 14, 2018 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ChrisRackauckas commented Feb 14, 2018

Choose a reason for hiding this comment

onoderat Feb 14, 2018 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

codecov bot commented Feb 14, 2018 • edited

Codecov Report

ChrisRackauckas commented Feb 14, 2018

ChrisRackauckas commented Feb 14, 2018

ChrisRackauckas commented Feb 14, 2018

onoderat commented Feb 14, 2018

ChrisRackauckas commented Feb 14, 2018

onoderat commented Feb 14, 2018

ChrisRackauckas commented Feb 14, 2018

onoderat commented Feb 14, 2018

ChrisRackauckas commented Feb 14, 2018

ChrisRackauckas commented Feb 14, 2018

onoderat commented Feb 14, 2018 • edited

ChrisRackauckas commented Feb 14, 2018

onoderat commented Feb 14, 2018 •

edited

onoderat Feb 14, 2018 •

edited

onoderat Feb 14, 2018 •

edited

codecov bot commented Feb 14, 2018 •

edited

onoderat commented Feb 14, 2018 •

edited