PoC Heun's method

[mtsokol]

Hi!

I'm a master degree student and I've recently started learning SDE topics, and I came across your project.

For the purpose of learning I've decided to try contributing - in known issues it states that Stratonovich methods are missing so I thought about trying that: Here's a draft of Euler counterpart - Heun's method.

If you would have a bit of time for a code review and guidance and it won't be too difficult for me I would like to try contributing.

It's just a draft with example run in demo notebook so without tests yet. What do you think about it?

Do you already have a concrete plan for Stratonovich module? Or do you suggest other good-first-issue to start?

Thank you for any help.

[googlebot]

Thanks for your pull request. It looks like this may be your first contribution to a Google open source project (if not, look below for help). Before we can look at your pull request, you'll need to sign a Contributor License Agreement (CLA).

:memo: Please visit https://cla.developers.google.com/ to sign.

Once you've signed (or fixed any issues), please reply here with @googlebot I signed it! and we'll verify it.

---

What to do if you already signed the CLA

Individual signers

• It's possible we don't have your GitHub username or you're using a different email address on your commit. Check your existing CLA data and verify that your email is set on your git commits.

Corporate signers

• Your company has a Point of Contact who decides which employees are authorized to participate. Ask your POC to be added to the group of authorized contributors. If you don't know who your Point of Contact is, direct the Google project maintainer to go/cla#troubleshoot (Public version).
• The email used to register you as an authorized contributor must be the email used for the Git commit. Check your existing CLA data and verify that your email is set on your git commits.
• The email used to register you as an authorized contributor must also be attached to your GitHub account.

ℹ️ Googlers: Go here for more info.

[lxuechen]

Hi,

Thanks for the PR! I think this is a fantastic first step in getting the Stratonovich solvers in place!

Before we move forward, there are two things I want to mention: 1) In order for this code to get merged, you would need to sign a Contributor License Agreement (CLA) per the contributing to google repo policy as hinted by the bot. IIRC, there are two versions: one for individuals and one for corporates. Going through the procedure shouldn't take too long.

2) Even though this is a simple solver, I'd still recommend writing up some numerical tests to see it's convergence order. I think a reasonable test is that we find a particular SDE, and test the Stratonovich solver's solution with varying step size against the solution from an Ito solver with very small step size that's already implemented in the code base.

With a proper correction term, a Stratonovich SDE can be converted to Ito form, and hence solved by an Ito solver. We could treat the solution found by the Ito solver as the "ground-truth" and thereafter inspect the convergence order of Heun's method by plotting the errors vs step size on a log-log scale. Here's an example I wrote for the SRK solver.

We don't have to include the code for this particular numerical test, and it might even seem like a slog to do the testing. But I think this is definitely important in the long run, since the last thing any user might hope for is a bug in one of the solvers.

If you're happy to complete these two steps, I'd be more than excited to move forward with this PR. I will also leave some comments on the code when I start working at some point over the weekend.

Let me know if there's more that I could help, and thanks again for starting this line of work.

[mtsokol]

@googlebot I signed it!

[googlebot]

CLAs look good, thanks!

ℹ️ Googlers: Go here for more info.

[mtsokol]

@lxuechen Thank you for your answer.

That sounds great! Let me spend some time on diving into the source code - I'll prepare the test and push it. Should it be done as mentioned diagnostic example or as additional method in TestSdeint? Or both?

[lxuechen]

I think it'd be nice if we could see a plot on log-log scale just to inspect the rate. We don't necessarily need to include the code for the test, since I'm not sure how much extra stuff that test might need. If the test code isn't too long and roughly follows the format of the other source files in the diagnostics folder, I'd probably name it stratonovich_diagonal.py and place it in the same folder.

[lxuechen]

Here's an example of what a rate-inspection plot would look like. I did this for the diagonal Ito methods, and you can actually run the code in the diagnostics folder to obtain this for yourself. The numbers in the parentheses are slopes and correspond to the strong orders of different methods.

For Heun's method, you should expect to get something that's close to 0.5.

[mtsokol]

@lxuechen Hi! I've managed to push first version of the stratonovich_diagonal.py diagnostics file. I used problem from skr_diagonal.py and run/plot it against euler and analytical solution:

...and the log-log plot:

Shapes of heun solutions look good but they seem to be slightly shifted in few cases (Is it OK?).
Also slope is close to 0.5.

I would be thankful for feedback - I hope that implementation is correct: I looked at formula and existing implementations:

• StochasticDiffEq.jl
• sdeint

This is almost the same, except for dW - it above implementations they use the same dWn for both stages where here I needed to compute next I_k for the second stage to get these solutions. I'm not sure if hadn't got it wrong while reading those implementations.

I can also polish PR in terms of naming and structure.

Also I think about minor restructuring - Could SDEIto be removed in favor of BaseSDE as it only gives sde_type (ForwardSDEIto and ForwardSDEStratonivich will inherit from BaseSDE and provide sde_type there)? Then all problems from problems.py can inherit from BaseSDE and be solved with both methods by providing concrete sde_type when instantiated. WDYT?

[lxuechen]

Hi,

Thanks for the updates! I just checked the code and the solver seems good.

Now regarding the rate inspection, there's a technical issue which is for an Ito SDE and Stratonovich SDE, their analytical solutions would be different, even when their drifts and diffusions are the same. So I am guessing that's the reason to the rather uniform vertical difference between Euler and Heun as you vary the step size.

One way to properly do the error analysis is to use a conversion rule to convert the Strat SDE to Ito form and then use the Euler solver based on Ito's form.

Now that I think about it, it might require more work than I originally thought, and things can become messy. I'm currently working on setting up the Stratonovich utils and have plans of doing the rate inspection properly. On the other hand, I don't want to leave your contributions unnoticed. So here's a plan that I propose:

It would be nice if you could only include the following files

• torchsde/__init__.py
• torchsde/_core/base_sde.py
• torchsde/_core/methods_strat/__init__.py
• torchsde/_core/methods_strat/diagonal/__init__.py
• torchsde/_core/methods_strat/diagonal/heun.py
• torchsde/_core/sdeint.py
• torchsde/_core/settings.py

I can take it up from there and set up all the stuff that's relevant for the analysis.

What do you think about this plan?

[mtsokol]

@lxuechen Thank you for your feedback!

That sounds good to me - I deleted diagnostic file leaving only 7 files requested in PR.

So what can I work on next? Maybe I can work on next method - Milstein method for Stratonovich problems? (I will finish gdg_prod in ForwardSDEStratonovich for that purpose). If not this maybe something different?

[patrick-kidger]

Closing in favour of #24