Switch to LevyArea.jl

[frankschae]

See #458
https://github.com/stochastics-uni-luebeck/LevyArea.jl
https://arxiv.org/abs/2201.08424

To do:

• Need to check AD compatibility [ see https://github.com/Can the Levy area computations be made compatible with AD? stochastics-uni-luebeck/LevyArea.jl#6 ]
• Update defaults [defaults to using LevyArea, allows for selecting specific routine, still contains separate implementations for diagonal and commutative noise which probably should be handled in a follow-up PR.]
• Update tests [no path-wise convergence tests]
• Update citations [ docs should be updated once this is merged.]

[fkastner]

A note on the terminology: Wiktorsson's algorithm is only one algorithm to approximate iterated integrals (or their associated Levy area). Another one is the one from Milstein (often attributed to Kloeden, Platen, Wright) or the newer Mrongowius-Rößler algorithm.

If the diffusion coefficient of your problem satisfies a special commutativity condition (which includes for example additive, linear and diagonal noise), then the Levy area part cancels out (because of its antisymmetry). In this case we can replace the iterated integrals by their symmetric part 1/2WW^\top and the SDE schemes get a lot easier since we don't have to approximate the Levy area anymore. But this has nothing to do with 'Wiktorsson approximation'.

Also I'm not really sure why you have different routines for diagonal and commutative noise? They should be the same, but there are some slight differences in your implementation I don't understand.

[frankschae]

Yup, agreed. The names go back to the first attempt where the goal was to implement Wiktorsson's algorithm (plus handling the special cases efficiently where the Levy area approximation is not required). I already updated the one that's user-facing
IIWiktorsson -> IILevyArea (open for better suggestions and suggestions for the other ones).

Also I'm not really sure why you have different routines for diagonal and commutative noise? They should be the same, but there are some slight differences in your implementation I don't understand.

This is because we represent the diffusion matrix for diagonal noise by a vector, see
https://diffeq.sciml.ai/dev/tutorials/sde_example/

[fkastner]

This is because we represent the diffusion matrix for diagonal noise by a vector, see
diffeq.sciml.ai/dev/tutorials/sde_example

Wouldn't it be easier to just wrap the return value of g in a Diagonal? That way the intention is clearer and the solvers don't need to be special cased.
But ok, that was probably discussed when the interface was created and is set in stone now.

[ChrisRackauckas]

It comes from the methods specialized on diagonal only, which came first. We can probably cut down on code by wrapping in a Diagonal, but it doesn't make a major change to the utility of the package so no one has done the work. Someone probably should refactor that though.

[ChrisRackauckas]

For strong convergence (and adaptivity) we would have to come up with a scheme so that we could easily pullback the Brownian paths. In the diagonal and commutative cases this is easy because it is done with the Brownian bridge, and the diagonal 1.5 case can be handled with a psuedo-Brownian process because the iterated integral approximation turns out to just be essentially two Brownian processes. I don't know what would have to happen with the non-commutative case though since it takes n random numbers?

[frankschae]

The test failures (StochasticDelayDiffEq/SROCKC2) look unrelated to me.

[fkastner]

For strong convergence (and adaptivity) we would have to come up with a scheme so that we could easily pullback the Brownian paths.

What exactly do you mean by pullback here?

Re adaptivity: I'm fairly confident to say that you can't use the Levy area in an adaptive scheme that requires the rejection of steps. We thought a bit about that, but it would require some new ideas or a completely new scheme for the Levy area approximation.

[fkastner]

Is the distinction between in-place and out-of-place still necessary for the iterated integrals? Especially since LevyArea.jl currently doesn't provide that distinction either. There is no reason why the in-place-ness of the iterated integral computations has to be coupled to the in-place-ness of the problem itself, or am I missing something here?

[frankschae]

I would say we could either remove the inplaceness completely, i.e., by removing it also from all SDE solvers

StochasticDiffEq.jl/src/perform_step/low_order.jl

Line 449 in c802f6e

get_iterated_I!(dt, dW, W.dZ, Wik, integrator.alg.p, integrator.alg.c, alg_order(integrator.alg))

(For problems with mutating functions, the iip formulations were originally used to avoid allocations -- this doesn't work for the methods from LevyArea.jl.)
Alternatively, we can use this wrapper that offers no computational advantage but allows us to remove allocations for the other cases (IICommutative()) in principle.

[ChrisRackauckas]

What the reason for not having a preallocated form in LevyArea.jl? I presume the allocation here would impact at least multithreaded performance.

[ChrisRackauckas]

The SDDE test is just noisy. I need to fix that when I get back.

[ChrisRackauckas]

Continuing the discussion on fully inplace in stochastics-uni-luebeck/LevyArea.jl#7 . While it could hit performance, for now we can probably just skip it.

[fkastner]

Alternatively, we can use this wrapper that offers no computational advantage but allows us to remove allocations for the other cases (IICommutative()) in principle.

Yes, that was my suggestion. Basically remove the out-of-place functions. It should be possible to use in-place Levy area approximations (no matter how much of it is in-place) even for out-of-place problems, shouldn't it?

[ChrisRackauckas]

out of place would still be required for some forms of automatic differentiation and for use cases with immutable arrays and static arrays. So that still would have a place in the grant scheme of things.