Euler-Heun method

[mtsokol]

Hi!
Just euler-heun method.
Tried to run it but it hangs on BrownianInterval so related to #36

[mtsokol]

@patrick-kidger Hi! Just rebased - now it all works:

[patrick-kidger]

Other than the minor comments, LGTM. Thanks for implementing.
The fix for BrownianInterval is now in dev so you should be able to run the diagnostics. What do they look like?
EDIT: looks like we crossed paths. What do the diagnostics look like for different noise types?

[mtsokol]

@patrick-kidger I've rewritten stratonovich diagnostics to scalar noise type (found one issue). Should I commit those files also? (it's basically copy paste of diagonal except for problem) Or only post plots here?
(and the same question regarding additive noise type)

[patrick-kidger]

As you've got them, let's include the files as well. We'll probably be interested in running those later too.
(In fact the diagnostics all share a lot of common code that should be factored out somewhere, but that's not super important now.)

[mtsokol]

Ok, I included stratonovich diagnostics for additive and scalar.
For both basic plots look OK.

Rates plots look like this:

• scalar
  
• additive [EDIT](changed back to Ex3)
  

There's an issue with milstein solver - for scalar noise dimensions does not conform to compute Y_prime as g it's one additional dimension bigger (included fix but it's a quick one just to make it pass).

And looking at rates milstein for additive it's wrong but haven't figured it out yet. Btw both grad-free and default use gdg_prod_default - according to comment grad_free path isn't ready yet, right?

[lxuechen]

So for additive noise SDEs, the diffusion doesn't depend on the state. So the Ito interpretation is the same as the Strat interpretation, and the Milstein schemes should all have the same strong order as the Euler scheme.

[mtsokol]

Thanks for review - now I got it. Updated plots above - now it all looks right.

[patrick-kidger]

The above if statement is superfluous. Also, the strong order should be 0.5 for general noise.

[mtsokol]

Ah, sorry for that, fixed.
Also I've started working on a benchmark from [19] from the paper ("Smoking adjoints"). For now I've tried to follow exactly their formulas and they use Euler solver in log form - I need to compute g there but for adjoint_sdeint it fails as it's not supported due to efficiency. What can I do about it? Can g somehow be computed for adjoint_sdeint? (I've also tried to make it as a field in sde and access it from _base_sde from solver but as in forward it's wrapped: ForwardSDE(base_sde) in backward I think it's AdjointSDE(ForwardSDE(base_sde)))

[patrick-kidger]

So back when that paper was written, "adjoint" meant "backprop through the solver".
These days, confusingly, we've appropriated the term to mean "don't backprop through the solver, and solve a differential equation instead".
As a result that paper is basically just about the now well-understood advantages of using backpropgation over forward-mode differentiation. (At least in this context.)

So they don't actually use the adjoint SDE at all, and what they do is equivalent to just using our sdeint in the normal way.

[mtsokol]

Yes, I got it while reading the paper. What I'm trying is to extend that benchmark from "smoking adjoints" with stochastic adjoint from here so I want to do it forward vs sdeint vs adjoint_sdeint.
And due to the fact that in custom log_euler solver that I created I'm using g the described error occurs.

[patrick-kidger]

I'm not sure what your distinction between forward and sdeint is.

You shouldn't need to create a custom solver I don't think - IIRC in that paper they use an SDE to model log-rates; solve that with Euler--Maruyama; take an exponential if you want the rates themselves.

Likewise, you shouldn't need to use g in the adjoint: the standard EM scheme only needs g_prod.

Does that make sense / am I understanding you correctly? If not perhaps can you provide some references to where in the SA paper you're looking? (There's a couple copies of that paper out there, I tend to look at this one - https://www0.gsb.columbia.edu/faculty/pglasserman/Other/RiskJan2006.pdf)

[mtsokol]

Yeah I was reading exactly that one.

By forward method I mean computing dV/dS with (sdeint(S+h) - sdeint(S-h)) / 2h, by sdeint ("adjoint" in "smoking adjoints", "backprop through the solver" here 😄) dV/dS is sol = sdeint(S); torch.grad(sol, S) and last one adjoint_sdeint - dV/dS is sol = adjoint_sdeint(S); torch.grad(sol, S). (I hope I didn't mess it up). Also I've got PoC and they yield same gradient result.

Code is similar to this (as forward vs sdeint): https://github.com/mgroncki/AdjointMonteCarlo/blob/master/PlainVanillaOption/blackMC.cpp

Thanks for help! I will try it that way. To explain myself the log-Euler solver was used in paper and also in references in Monte Carlo Methods In Financial Engineering (page 398) directly for that model.

[mtsokol]

But I think I need a bit more time to fully figure all these methods out on paper.
The chapter 7. of that book also makes comparison of two methods: the finite-difference approach defined as (Y(S+h) - Y(S-h)) / 2h and pathwise mehtod (which is the one compared in "smoking adjoints").

[patrick-kidger]

This LGTM. @lxuechen happy for me to pull the squash+merge trigger?

[lxuechen]

The code in general is okay. I'd recommend using a linter with the character limits and blank line between imports/functions enforced.

Thanks again for the putting in the effort in implementing these methods.

[mtsokol]

Thanks for review! My main PL is Scala and python integration in my VS code barely works but I will setup that linter options to avoid those comments in the future.

[mtsokol]

Just one moment - milstein naming should be fixed also in diagonal and scalar files.

[mtsokol]

@lxuechen OK, now it's ready.

[mtsokol]

P.S. (asking out of curiosity) Do you have a roadmap for merging dev into master and cutting new release?

[patrick-kidger]

@mtsokol
Loosely speaking yes, although @lxuechen and I haven't fixed things in detail.

I think we'd like to tick off the remaining things in #15. Fortunately, I think we've now done all the complicated ones, and it's mostly just grunt work remaining.
My hope is for a new release this month. (I'm only speaking personally there; Chen and I haven't discussed release schedules yet.)

[lxuechen]

@mtsokol @patrick-kidger

I agree mostly with Patrick. Aside from the technical problems in #15, on a high level, I think the most crucial missing things are

• documentation (docstring, README, additional doc site generated from docstring?)
• tests, tests, and tests
• benchmarks
• continuous integration
• backward comp fixes

I think we should be able to pull off another release in this month, if things go well. Though, I also think we shouldn't rush things by too much and only merge when we're genuinely ready.