Milstein scheme for SME and a cool example for testing of a stochastic solver

[Vutshi]

Hi all,

I'm not sure how to do pull request in this case so let it be a new Issue.

I have implemented Milstein method for solving SME for homodyne detection scheme. For my problems it works much better than the simple Euler-Maruyama method. In case of a single Wiener increment all what is required is the new rhs function (see notebook attached). For multiple Wiener increments the provided infrastructure of QuTiP is not enough. One needs to provide all of the A_ops to the rhs. The modified stochastic.py file is also attached (from qutip.stochastic import smesolve_mil).

The attached Notebook contains an example of SME with an analytical solution. Which allows to test smesolver for a single quantum trajectory.

The notebook and the modified stochastic.py is here:
http://db.tt/dJNjGo3g

P.S. Due to nonlinearity in the stochastic term of quantum SME the simple Euler-Maruyama scheme doesn't work properly sometimes. Even Milstein struggles. From my experience semi implicit methods help. But in order to implement it in a convenient way one needs to modify QuTiP approach to the problem. Maybe instead of rhs function it would be better to have a function which returns \rho_{n+1} instead of d\rho. It could be more flexible.

[jrjohansson]

Hi Vutshi

Thanks a lot for this contribution! It's looking very interesting indeed. I'll be looking into the details of your code, notebook and your paper in the coming days, but from a first look it's looking great! I've been working a lot on the stochastic solvers the last few months, but it is still very much under development. The milstein method has been on my todo list for a long time, but I've not gotten around to work on it yet, so your contribution is very timely and most welcome!

I have tried to design the API in the most general way I could think of while still keeping the complexity from growing too much. Defining an SME with its RHS function, D1 and D2 functions has worked for most problems I considered so far, but I'd be happy to discuss further if this could be changed to something more flexible, like a function for calculating rho_{n+1} like you suggest. We hope to have the stochastic solver API ready and mature enough to include in the next qutip release sometime this fall, so it would be good to have these things tested and thoroughly discussed before then. But until the stochastic module is officially in the QuTiP API we still have a lot of freedom to shake things up if necessary.

Also, to keep complexity and code duplication down I wanted to follow the pattern

smesolve (choose solver by argument) -> smesolve_generic -> specific method and SDE by rhs, d1, d2 selection

So it would be good if smesolve_mil and smesolve_impl in your code could be combined with smesolve, and if smesolve_implicit and smesolve_milstein could be combined with smesolve_generic, etc. That would reduce a lot of code duplication I think. Of course, if smesolve_generic needs to be changed to accommodate this that would be totally fine with me.

It would be great if you could do a pull request of this to the qutip master, then we could discuss details in the PR. You'd need to rebase your branch to the current master tip though, because things already changed a bit since the version you based your additions on. If you have troubles with this or no not know how to do it I could also try to merge your modifications into the qutip master, it's totally up to you.

[Vutshi]

What I did in stochastic.py is a totally redundant mess. This is the reason I didn't do the pull request. I feel that an elegant way to add the Milstein with multiple Wiener increments (as well as any other more advanced methods) is to switch to a RHS function which is provided with all A_ops and which returns rho_{n+1}. Since it involves changing of API I want to discuss it with you first.

[jrjohansson]

Hi again

I've looked at your notebook and code in some more detail now, and I have attempted to make the change in how the rhs function behaves. I had not seen an example where the contributions form different stochastic collapse operators depended on each other before, but I now see the need for having this possibility from looking at your implementation of the Milstein implementation with multiple Wiener processes.

I have created a new branch in qutip called sme-rhs-restructuring (see commit 0b4a1a7 and https://github.com/qutip/qutip/tree/sme-rhs-restructuring) where I have experimented with rhs functions that internally loop through the A_ops list. In simple cases it only loops through the list and calculates a contribution for each stochastic collapse operator, but now it is also possible to do something more sophisticated, like in your milstein solver. I've edited your notebook so that it uses the standard (now modified) smesolve and your rhs functions for the milstein methods.

http://nbviewer.ipython.org/6153688

Note that I still only use the rhs function for the contributions due to the stochastic collapse operators (both the deterministic and stochastic part), and it does not include the unitary part or the deterministic collapse operators (if those are given as well), so it is not quite what you suggested when proposing that the rhs should calculate \rho_{n+1}. Do you think this approach is reasonable, or do you see any problems with it?

Regarding your milstein implementation, it seems that it is only valid for homodyne detection, since your have hardcoded the d1 and d2 function into the rhs function. I've tried to separate the definition of the SDE (d1 and d2) from the implementation of the SDE integrator (rhs), and naively it seems like you could use the d1 and d2 functions in rhs_milstein to make it possible to use it with hetrodyne detection as well. Is there any difficulties in generalizing the milstein method and parameterize it in terms of d1 and d2 functions in this way?

I've also tried to use the smesolve_imp and smesolve_milstein with explicit=False in your original notebook, but I do not get resuls that agree with the other methods if I use these. Looking at the code for smesolve_imp, you do an inversion of some variant of the liouvillian to calculate Lsparse, what is the purpose of that calculation? And finally, I see you use expm in _smesolve_single_milstein_implicit: I think this should be avoided because it will be very computationally demanding for larger systems.

If you are satisfied with how the new rhs functions work I'll merge the branch sme-rhs-restructuring into qutip master and then feel free to submit a pull request with your milstein rhs implementation.

[Vutshi]

Hi,

I think your new implementation of API should be flexible enough to incorporate any explicit scheme. The reason I think it is useful to calculate \rho_{n+1} is to apply a semi implicit method. This is what is done by Lsparse. I tried two ways of doing it. First, one can calculate (1 - (L+d1)_dt)^(-1) and apply it to to the stochastic part of the equation. Another way is to calculate the deterministic evolution expm[(L+d1)_dt] and also apply it to the stochastic part. Both of these methods work for linear deterministic part like in the homodyne detection. I've got good results with these semi-implicit methods on a problem with a few qubits. Due to a big nonlinearity for some parameters regime it really helps to apply the true deterministic evolution (or its approximation). For another regime it is better to use explicit method. In principle it should be good to adjust degree of implicitness. Anyway I'm still looking into the problem and comparing different possibilities. I hope to produce a new notebook with tests soon.

This Milstein realisation is compatible in principle with the heterodyne detection. I'm not sure it is possible to express everything in terms of d1 and d2 only. The point is that I need to use derivatives of d2. I don't really see a problem to have separate RHSs for homodyne, heterodyne and photocurrent detections. I'll think a bit more about unification of homodyne and heterodyne and then I'll submit a pull request.

Regarding expm for large systems. We need to calculate it only once, it works for small systems so I think it is a good option to have.

[jrjohansson]

Hi

OK, I see. For implicit SDE solvers there is and advantage in having the rhs function calculating the state at the next time-step. I hesitated doing it like that, because I didn't want to mix up the definition of the stochastic SDE contribution with the deterministic part, which sometime in the future should use qutip's existing framework for time-dependent hamiltonian and collapse operators. Having the stochastic rhs function exposed to all that is a bit unnecessary and makes things more complicated than it needs to be. However, I think I found a reasonable compromise in that I created new _rhs_psi_deterministic and _rhs_rho_deterministic where the current deterministic contributions are calculated. These functions can now be called from a stochastic rhs implementation, like this

dpsi_t = _rhs_psi_deterministic(H, psi_t, t, dt, args)

In this way the implementation of a stochastic rhs function is relatively isolated from how the deterministic part is calculated (which right now is trivial, but which could become more complex when time-dependent systems are implemented). However, it also gives all the flexibility needed in the rhs function, and these deterministic rhs functions do not need to be used if it is not suitable in a particular stochastic rhs implementation.

I've added the changes to the rhs behavior to the https://github.com/qutip/qutip/tree/sme-rhs-restructuring and updated the gist with the variant of your notebook

http://nbviewer.ipython.org/6153688

Note that a few other things, like the rhs function signature, had to change in the process. I hope that this updated API will be flexible enough to let you implement implicit solver.

Regards rhs implementation and d1,d2 function: Yes, I understand that not all possible rhs schemes might be compatible with the parameterization with d1 and d2 functions, but when it is possible it has the great advantage that the user do not need to worry about how the rhs is implemented, only define the SDE in terms of d1 and d2, and then possibly select a rhs solver using the solver argument to smesolve. This will not work in general for solvers like the milstein scheme, since it requires an analytical derivative (unless it can be evaluated numerically?), and in those cases it would be sufficient to implement problem specific rhs functions (like rhs_rho_milstein_homodyne etc).

Regarding the use of expm: OK, I agree that it could be a nice method to have to be used on smallish system. However, since we are stepping with a small timestep dt, shouldn't it be sufficient to use an expansion of of expm as usual? Perhaps a second-order expansion would do if the first order isn't sufficient. Doing a full expm is fine too I guess, but if it can be avoided then the same method could be used on larger systems too.

Looking forward to seeing your new notebooks and the PR!

[Vutshi]

Hi Robert,

I was thinking about heterodyne and homodyne detection unification. It seems to me that the way how the heterodyne is done now in QuTip is a bit overcomplicated. Why don't we just produce two jump operators $s/sqrt2$ and $-is/sqrt2$ out of one heterodyne $s$ and run the homodyne scheme with these two operators. The Milstein will work automatically. Do I miss something?

[jrjohansson]

Hi Vutshi

Yes, it is perfectly possible to implement the heterodyne detection using two homodyne detections with 50% efficiency, like you suggest. For the Milstein solver it might be the best, or at least the quickest, way to implement heterodyne, but I don't really agree that the current method is overcomplicated. It is in fact basically the same as the method you suggest, so I don't see that either would be more complicated or simpler than the other. The only issue is whether the two stochastic increments are divided up in several stochastic collapse operators or if the d2 function internally takes care of the both stochastic increments for a given stochastic collapse operator. The complexity is just shifted from the d2 function to somewhere else, which in general doesn't simplify anything. However, since the milstein solver you submitted is written so that it only support one increment per collapse operators, then there might be a real advantage of splitting the heterodyne process into two homodyne processes.

I've tried to make the qutip stochastic solver API as general as I could (and it is still a work in process), so that it will be as flexible as possible for implementing custom types of SMEs. Having support for multiple stochastic increments per collapse operators seems to be useful in certain applications. Although such SMEs could probably always be rewritten as multiple collapse operators with single stochastic increments, it might not always be the most natural way to define the SME. The heterodyne detection is one example of this, and it can be formulated in both ways, but I want the qutip SME API to work with both methods (at least with the basic euler solver, not necessarily with every solver we implement). However, that doesn't mean that the implementation of heterodyne for a particular solver has to use one way or the other. We should just document which solvers support multiple increments per operators and which solvers don't.

So if it is easier to get the milstein solver working for heterodyne by simulating two homodyne detections, then let's go ahead use that method for that solver. The required preprocessing can be done in the smesolve function, or something like that.

[jrjohansson]

Resolved by PR #41.