Simulating large nontrivial SDE systems: algorithms, reliability and speed

Specific Domains Modelling & Simulations
, ,
1

In short, I am interested in the simulation of large systems of nontrivial stochastic differential equations as they arise when applying phase-space methods to interacting quantum systems.

As an intro example: when one uses the so-called “positive P representation” on the two-site Bose-Hubbard model, this would result in eight nonlinear equations for the deterministic part, along with four independent Wiener processes that then enter the system as multiplicative noises, i.e noise_rate_prototype=zeros(8, 4) .

Note: I am not asking how to solve this specific system.

Rather, I would like to have some input on the more general questions:

  1. What algorithms are suggested for systems that have both, nonlinearity and multiplicative noise? (For the two-site Bose-Hubbard, I have tested ImplicitEM(), with satisfactory results in ‘unproblematic’ parameter regimes.)

  2. Numerically, how large can I make such systems and still get a reliable solution in finite time? For the combination of Bose-Hubbard and positive P, for instance, the degrees of freedom will have a factor of four w.r.t. the number of sites, and I will have half as many independent Wiener processes. Say I have 128 sites, that will make 512 degrees of freedom and 256 Wieners, combined in a coupled system of equations of the type described above. My naive intuition suggests that, among other things, one will have to create ever larger ensembles when the number of noises goes up. The nonlinearity will add further complications, obviously. Hence the issue of speed vs. accuracy seems to be of great concern here, too.

Thanks in advance for any help!

3 Likes
2

Great question with lots of answers to explore.

General Overview

The current state of where we are at is well-summarized in this preprint:

Not all of the algorithms are in there, because for example you’ll notice that there are some algorithms here which have no source:

http://docs.juliadiffeq.org/latest/solvers/sde_solve.html#Stiff-Methods-1

I will get those written up, but they are adaptive and for the large noise case since methods like ImplicitEM() are still largely unstable when the noise term is sufficient large (since implictness on SDE solvers is semi-implicit only in the drift term due to the non-invertability of standard normal random variables, though there are some truncation schemes that can be done and this is quite a rabbit hole I would like to go down but not right now).

But okay, so the state is basically that if you have an additive noise SDE there’s a ton of stuff we can do, and that paper shows how to convert an affine noise (or multiplicative noise) SDE into an additive noise SDE (under certain conditions for invertability of the noise term). If it’s a stiff diagonal noise SDE, in many cases SOSRI is doing better than the (semi) implicit algorithms which is showing that there may be a deficiency in their design, something that has been noted before:

We have plans (and students!) to try out alternatives like the SROCK methods, track it here:

What it means for you

Well there’s a few more things to mention. When you say four independent Wiener processes, do you mean real-valued Wiener processes or complex-valued Wiener processes? If they are complex valued, then your SDE isn’t diagonal but actually commutative (since multiplication mixes the real and complex valued terms) and so the diagonal noise methods will have order loss to 1.0. This is another reason to say that by far the best technique will be to convert to an additive system using the Lamperti transform mentioned in the paper, and then applying SOSRA or SKenCarp.

One thing we do need to add to our implementation though is the ability to tell it how many random numbers you want. I haven’t seen this factor into profiles though, and the diagonal/additive noise methods are so much more efficient than the general SDE methods like ImplicitEM that I would be heavily surprised if using noise_rate_prototype=zeros(8, 4) is actually helpful here. On the topic of profiling…

Please donate an example problem

Most of the testing and development has centered around a set of 20 highly stiff biological SDEs because someone gave this nice problem to me:

https://github.com/JuliaDiffEq/DiffEqProblemLibrary.jl/blob/master/src/sde_premade_problems.jl#L172-L318

Please feel free to donate an implementation of this problem for us to benchmark and profile on. That would be very helpful for algorithm development since for SDEs it’s not really known what all of the properties of equations that are important for real-world performance, so testing directly on problems related to the application can be a good idea.

Also, if you cannot convert it to an additive noise problem due to invertability issues and it’s a complex Wiener process, that would be a great real-world commutative noise example and I would really like to have an implementation of that for testing some algorithms in development!

9 Likes
3

Thank you for the detailed answer! I will post the problem here after I have digested your input!

4

Here is the example implementation of the problem. When I said “four independent Wiener processes”, I actually did mean real increments, although I realize now that I probably should not have called the noise “multiplicative”… It is, in some sense, as you can see below, but seemingly not in the convenient “diagonal” form that would lend itself immediately to a Lamperti transform, as far as I can see. This was probably expected since simulating interacting quantum systems should not be that easy.

Concerning the parameters I have set below, it should be noted that they are very small: one certainly wants the solution of these equations for either large U or N, and possibly also for cases where both are not small. Furthermore, the deterministic system has an instability in terms of the dimensionless parameter \Lambda = NU/2J, as well as a so-called “self-trapped” regime for large initial values of z. In this sense, the simple two-site model is already very rich and a good playground to start off from.

As I said, I already have obtained some results for the parameters below with ImplicitEM(), but I certainly did have problems stretching that approach to the more interesting parameter regimes…

# some arbitrary evolution time
t_max = 1.

# hopping
J = -2*pi*1.0/2.

# interaction (this is a very weak value, much larger ones definitely desirable)
U = 2*pi*0.005

# (very small) number of particles (should be approx. conserved)
N = 10.0

# deterministic part (linear terms prop. to J, nonlinear terms prop. to U)
function sys(du,u,p,t)
  du[1] = J*u[6] + U*u[1]^2*u[7] + 2*U*u[1]*u[3]*u[5] - U*u[5]^2*u[7]
  du[2] = J*u[5] + U*u[2]^2*u[8] + 2*U*u[2]*u[4]*u[6] - U*u[6]^2*u[8]
  du[3] = -J*u[8] - 2*U*u[1]*u[3]*u[7] - U*u[3]^2*u[5] + U*u[5]*u[7]^2
  du[4] = -J*u[7] - 2*U*u[2]*u[4]*u[8] - U*u[4]^2*u[6] + U*u[6]*u[8]^2
  du[5] = -J*u[2] - U*u[1]^2*u[3] + 2*U*u[1]*u[5]*u[7] + U*u[3]*u[5]^2
  du[6] = -J*u[1] - U*u[2]^2*u[4] + 2*U*u[2]*u[6]*u[8] + U*u[4]*u[6]^2
  du[7] = J*u[4] + U*u[1]*u[3]^2 - U*u[1]*u[7]^2 - 2*U*u[3]*u[5]*u[7]
  du[8] = J*u[3] + U*u[2]*u[4]^2 - U*u[2]*u[8]^2 - 2*U*u[4]*u[6]*u[8]
end

# noise strength (related to interaction)
sigma = sqrt(U)/2

# nondiagonal (?) multiplicative noise: the Wiener increments
# for du[1, 1] and du[5, 1] are the same, and similarly for
# the remaining three pairs (du[2, 2] and du[6, 2] etc.)
function σ_sys(du,u,p,t)
    du[1, 1] = -sigma*(u[1] + u[5])
    du[2, 2] = -sigma*(u[2] + u[6])
    du[3, 3] =  sigma*(u[3] - u[7])
    du[4, 4] =  sigma*(u[4] - u[8])
    du[5, 1] =  sigma*(u[1] - u[5])
    du[6, 2] =  sigma*(u[2] - u[6])
    du[7, 3] =  sigma*(u[3] + u[7])
    du[8, 4] =  sigma*(u[4] + u[8])
end

# so-called population imbalance z = (N_1 - N_2)/(N_1 + N_2)
z = 0.1

# initial population derived from z
N_1 = 0.5*N*(1+z)
N_2 = 0.5*N*(1-z)

# SDE problem where initial phases are set to zero, hence real initial values
prob_sde = SDEProblem(sys, σ_sys, [sqrt(N_1), sqrt(N_2), sqrt(N_1), sqrt(N_2), 0., 0., 0., 0.], (0.0, t_max), noise_rate_prototype=zeros(8, 4))
5

Let me remark a couple of other things.

Concerning the people who invented the mapping from quantum to classical description I am using above, they have provided an algorithm here. There is also a software package related to that called xmds. (I can only post two links for a start?) The algorithm is still in use in more recent work.

6

Robust Algorithms for Solving Stochastic Partial Differential Equations - ScienceDirect is fine. The algorithm they build in to 2.9 is an implicit form of RKMilCommute(interpretation=:Stratanovich). Does your model satisfy the commutivity relation? It wouldn’t be hard to implement the implicit form (it’s mostly pasting together ImplicitRKMil and RKMilCommute and adding tests) if you need it.

http://docs.juliadiffeq.org/latest/solvers/sde_solve.html#Special-Noise-Forms-1

If it does then that’s likely the method you want to be using. Or for large sigma you might want to use ISSEM because the semi-implicit methods (like ImplicitEM and the ImplicitRKMilCommute from the paper) are not stable for large noise eigenvalues, while step-splitting corrects for this. I think what you really want is an ISSRKMilCommute.

So please do check and if it does satisfy the commutivity equation I can add these algorithms since this is really what you need and I want to start researching better methods for this case anyways.

The rest of the paper is just the spatial discretization which is part of the SDE definition so I’ll leave that for you.

1 Like
7

Okay, so it was quite easy to see that the noise is indeed commutative. I get that it makes a lot of sense then to exploit this property algorithmically.

Thanks again for your great support!

8

Yup, it means the Wiktersson approximation of the iterated integrals decomposes and then is much simpler to compute.