0.2.0 TODO list

[patrick-kidger]

To-do list so far:

Brownian motions:

• Switch to interval interface.
  • Brownian
  • BrownianInterval
  • BrownianPath (both Python and C++)
  • BrownianTree (both Python and C++)
  • calls in adjoint
  • calls in solvers
  • calls in examples/demo.ipynb As of Dev kidger2 #19 support single-point evaluation
  • calls in tests/problems.py As of Dev kidger2 #19 support single-point evaluation
  • calls in tests/test_brownian_path.py
  • calls in tests/test_brownian_tree.py
• Levy area support:
  • BrownianInterval
  • BrownianPath (Brownian unification #61)
  • BrownianTree (Brownian unification #61)
  • Removed 'trapezoidal_approx' as option from SRK solvers (see Replace trapezoidal_approx with integrated Brownian bridge for speed #12)
  • Removed 'trapezoidal_approx' from examples (see Replace trapezoidal_approx with integrated Brownian bridge for speed #12) (Done in Dev kidger2 #19)
• Tidy up __init__:
  • BrowianPath (Chen) (Add support for Stratonovich adjoint #21)
  • BrownianTree (Brownian unification #61)
• Misc:
  • Handle t outside of t0, t1. (Patrick) (Global entropy doesn't determine sampled values outside of [t0, t1] for BrownianTree #9, Calling BrownianTree with t > t1 returns None #10, BrownianInterval tests + bugfixes #28)

Solvers:

• Not split up by noise type any more, to the extent that's reasonable.
• Add Midpoint solver
• Add log-ODE solver (Add log-ODE scheme and simplify typing. #43)
• Allow picking the solver for the adjoint pass.
• Remove integrate_logqp and interp_logqp
• Fix step sizes leaving the region of integration. (Patrick) (BrownianInterval tests + bugfixes #28)

SDEs:

• Add gdg_prod for the adjoint. (Chen)
• Add general-noise adjoint (Chen)
• Tidy up the directory structure for the adjoints a bit.
• Add checks in the adjoint for Ito/Stratonovich correction terms. I think the checks should depend on both the sde_type and the method.

Misc:

• Bumped version number, python_requires.
• diagnostics/ito_scalar.py appears to be using huge amount of memory, and is also using a diagonal (and thus not compatible) problem.
• Update documentation. (Added documentation #71)
• Better adaptivity checks. (e.g. general-noise Heun should be possible with adaptive stepping, I believe.) (Dev methods fixes #73)
• Support logqp (Dev logqp #75)
• Change sdeint, sdeint_adjoint to automatically select the probably-best solver for a given noise type / SDE type combination. (Rather than just always using SRK.) (Dev methods fixes #73)
• SDE-GAN example

Tests:

• Fix broken tests test_sdeint.py, test_adjoint.py, test_adjoint_logqp.py (Chen) (Add support for Stratonovich adjoint #21)
• Write tests for BrownianInterval (Patrick) (BrownianInterval tests + bugfixes #28)
• Add repetitions back to test_brownian_interval.py::test_normality_conditional (BInterval - U fix #44)
• Write distribution tests for space-time Levy area. (BInterval - U fix #44)
• Update test_sdeint.py to pytest. (Dev methods fixes #73)
• Add tests checking that every solver works. (Dev methods fixes #73)
• Fix flakiness+slowness of BInterval tests (Fixed BInterval test flakiness #76)

Bugs to fix:

• Correction a from space-time Levy area to Levy area is of the wrong shape. (Patrick) (BrownianInterval tests + bugfixes #28)
• Davie's approximation is currently missing its correction term from the Brownian parabola. (Patrick) (BrownianInterval tests + bugfixes #28)
• Use the correct noise term for I_k0 in SRK.
• Default values of method='srk' and levy_area_approximation='none' are incompatible with one another. (Done in Dev kidger2 #19)
• If bm isn't passed then it defaults to being a BrownianPath on the CPU, potentially causing device errors later. (Done in Dev kidger2 #19)

[patrick-kidger]

@lxuechen Benchmarks look rather encouraging for the BrownianInterval. e.g. see below. (Every benchmark exhibits the same behaviour.)

I would note that the number of time steps is very low in the benchmarks (only 100). I expect that larger timesteps will result in the BrownianPath overtaking the BrownianInterval.

[lxuechen]

This is encouraging! Would you have an idea as to where the gain is coming from?

I would guess it's coming from 1) using __slots__, and 2) explicitly modeling increments.

How much different are the plots when these bms are plugged in the solvers?

[patrick-kidger]

I can think of a couple other things for certain cases:

• For sequential access: the algorithm used in BrownianInterval is reasonably smart about where it starts searching; IIRC BrownianPath doesn't use any hinting about where to start searching.
• For random access with BrownianPath: not having to do list inserts, which are O(n) on average.
• For the GPU: not to having to do a .to() copy and instead initialising right on the GPU via device=. TBH this is really just a bug in the other implementations though.

The equivalent plots for the solvers look pretty much identical; I just picked one at random.

[lxuechen]

Several side comments:

Switch to interval interface.

This breaks backwards compatibility. I'd be nice if we could have tb as an optional argument for BPath and BTree for now.

Removed 'trapezoidal_approx' as option from SRK solvers

This line might be wrong if the second output of bm is H in James' thesis, since I_k0 isn't exactly the space-time Levy area.

A side comment is that I don't totally support merging the SRK methods for now. But I think we should be able to get rid of the tableuas folder if you find it disturbing.

Also the docstring needs to change for that file.

[patrick-kidger]

Several side comments:

Switch to interval interface.

This breaks backwards compatibility. I'd be nice if we could have tb as an optional argument for BPath and BTree for now.

Done.

Removed 'trapezoidal_approx' as option from SRK solvers

This line might be wrong if the second output of bm is H in James' thesis, since I_k0 isn't exactly the space-time Levy area.

I'll add that to the bug list!

A side comment is that I don't totally support merging the SRK methods for now. But I think we should be able to get rid of the tableuas folder if you find it disturbing.

The SRK methods actually aren't merged at the moment; for now at least I've taken the reasonably straightforward approach of unifying as much of the easy code as I can and leaving their steps separate. The tableuas folder doesn't bother me. :)

Also the docstring needs to change for that file.

Done!

[lxuechen]

I can think of a couple other things for certain cases:

• For sequential access: the algorithm used in BrownianInterval is reasonably smart about where it starts searching; IIRC BrownianPath doesn't use any hinting about where to start searching.
• For random access with BrownianPath: not having to do list inserts, which are O(n) on average.
• For the GPU: not to having to do a .to() copy and instead initialising right on the GPU via device=. TBH this is really just a bug in the other implementations though.

The equivalent plots for the solvers look pretty much identical; I just picked one at random.

I see. For GPUs the copying part makes sense. The list inserts argument doesn't seem to be the most convincing to me, since blist should have O(logn) complexity.

[lxuechen]

Overall seems like awesome progress!

[patrick-kidger]

Good point about the blist.

Thanks. I'm expecting to try and tackle most of the bugs first. I've just seen your email so I'll continue responding over there.

[lxuechen]

I've left the test and examples alone as I suspect you'll be much more familiar with how those should go than I am.

I'd recommend we first test BrownianInterval and check for its correctness. This actually should have been done before rewriting the solvers. The problem is that now, I couldn't run the rate inspect diagnostics, since the BInterval is raising errors.

I think all the 'SDE' things are the bits that I'm least familiar with, and that you're most familiar with. Can you handle them?

I can do this.

[patrick-kidger]

I've left the test and examples alone as I suspect you'll be much more familiar with how those should go than I am.

I'd recommend we first test BrownianInterval and check for its correctness. This actually should have been done before rewriting the solvers. The problem is that now, I couldn't run the rate inspect diagnostics, since the BInterval is raising errors.

One of the reasons I had it in a separate branch originally!
What error is being raised?

I think all the 'SDE' things are the bits that I'm least familiar with, and that you're most familiar with. Can you handle them?

I can do this.

Excellent, thankyou!

[patrick-kidger]

If my derivation is correct then I think the term needed in the SRK should be:

I_k0 = (t - s) * (H_{s,t} + 0.5 W_{s,t})

But it's also 1am where I am. Does that look about right to you?

[lxuechen]

If my derivation is correct then I think the term needed in the SRK should be:

I_k0 = (t - s) * (H_{s,t} + 0.5 W_{s,t})

But it's also 1am where I am. Does that look about right to you?

Yes, this is right. I_k0 is exactly the U_{s, t} term in the note I sent you, and the conversion is eq (4) in that note.

[lxuechen]

I think, for the moment, I wouldn't worry too much about rewriting/refactoring more of the existing code.

The thing that I'm worried more about is whether the new BrownianInterval is resilient enough to stand through numerical tests. The idea here is that I'm seeing subtle issues introduced by this new data structure and the various rewrites, and I don't want to be moving forward before we properly test BrownianInterval.

This also relates to the habit of mine where I'd like to run the test suit every time I make a minor change, so that I know exactly what could go wrong.

[lxuechen]

Here are tests that I think would benefit us when we go forward for BrownianInterval (and even for BPath and BTree when I add in the L.A. approxs)

• Shape and dtype tests for all cases -- either the L.A. is returned or not.
• The method to, i.e. to a new device or dtype should transfer/convert all existing tensors -- including those in the cache if possible.
• Determinism of the implementation -- a test like this -- but slightly more complicated. It should involve checking that the spawned increments and L.A. approxs are the same even when the cache is tiny and things are forgotten eventually.
• Normality tests and moment matching tests: This one could be tricky. Is there a way to set the Brownian motion value at a specific time like this so that we could test the conditional moments?

I think once these are in place, I'd have much more confidence in going forward.

[lxuechen]

Another note is that is there a way that we can simplify the interface, as right now the user is basically forced to manually supply levy_area_approx='davie'/'foster'/'space-time' just in order to use SRK. Using a high order scheme shouldn't be this difficult!

Also, does Davie approximation have any advantages over Foster? I thought Foster's version is more accurate, at the cost of one extra random number generation, which is totally reasonable considering the costs of other components.

[patrick-kidger]

Your concerns about tests are reasonable - bugs and tests are definitely the next things to handle.

I agree that having to specify levy area to use SRK should be fixed; this is already on the bug list. I'm not sure what the best way to handle this is.
I think my preferred fix would be to upgrade the default sdeint(..., bm=None) case to initialise a Brownian object with the appropriate levy area, and regard the explicit bm=something case as more of a power user's case, for which it's up to them to do things correctly.
I also like the ability to then switch between davie/foster approximation independently of the solver, in particular as better approximations may be introduced in the future.

The alternative option is slightly messier, but still doable -
Brownian* should be fixed to a particular choice of levy area approximation, as doing otherwise introduces headaches: what if it's queried with one choice, and then a different one?
Having the choice of levy area approximation be specifiable by the solver, rather than on __init__, then means that some sort of "deferred __init__" is necessary, which is called when Brownian* is passed into sdeint. I think that would only be used to set the levy area approximation flag, and nothing else, so this is definitely doable, but I'm not a fan of having partially initialised objects wandering around.

Care to make a judgement one way or the other?

Davie's approximation is slightly faster than Foster, by a bit more than just the random number generation (there's a bit of algebra too), but AFAIK there's not a big difference between them however you slice it. But for forward compatibility with even better approximations in the future, I'd avoid code layouts for which Foster gets baked in as the only option.

[lxuechen]

I think my preferred fix would be to upgrade the default sdeint(..., bm=None) case to initialise a Brownian object with the appropriate levy area, and regard the explicit bm=something case as more of a power user's case, for which it's up to them to do things correctly. I also like the ability to then switch between davie/foster approximation independently of the solver, in particular as better approximations may be introduced in the future.

I agree with this, and the plan seems nice.

Davie's approximation is slightly faster than Foster, by a bit more than just the random number generation (there's a bit of algebra too), but AFAIK there's not a big difference between them however you slice it. But for forward compatibility with even better approximations in the future, I'd avoid code layouts for which Foster gets baked in as the only option.

SGTM.

[patrick-kidger]

@lxuechen This is labelled as Ito-specific:

torchsde/torchsde/_core/base_sde.py

Line 53 in 133a745

class ForwardSDEIto(SDEIto):

Is there any part of that which actually relies on that though? It looks like it's just general drift/diffusion operations.

[lxuechen]

Is there any part of that which actually relies on that though? It looks like it's just general drift/diffusion operations.

This is removed in #20 and #21.