Steadystate function in the 3.0.0 version

[dan hu]

Hi everyone,

I installed the qutip of the version of the qutip-master (3.0.0 version) and found something very wired about the steadystate function. The attached code is aimed to calculate the steady-state g2 function for optoemchanical system in the strong coupling regime where g0=2wm. But it output a negative value of g2 =  -566.8944. Is that possible? or because the subspace (nc, nm)=(3, 160) is too small here?

The parameters  I used are:

nc = 3 # dimension of the cavity mode

nm = 160 # dimension of the mechanical mode

wm = 2*pi

g0 = 2.0*wm

kappa = 0.005*wm

gamma = kappa/3

nth=20 # mechanical temperature

cdr= 0.001*kappa

H = delta * Nc + wm * Nm - g0*Nc*(bm+bm.dag()) + cdr*(ac.dag() + ac)

with collapse operators

c_op_list = []

r1 = gamma * (nth + 1)

r2 = gamma * nth

c_op_list.append(sqrt(kappa)*ac)

c_op_list.append(sqrt(r1)*bm)

c_op_list.append(sqrt(r2)*bm.dag()) 

rhoss=steadystate(H,c_op_list)

g20=expect(ac.dag()*ac.dag()*ac*ac, rhoss)/expect(ac.dag()*ac, rhoss)**2

The code run about 15 mins. 

Thank you for your help and your time!

Best,

Dan

[dan hu]

By the way, it seems that it only happens for very small drivings amplitude cdr. If I change cdr = 0.01kappa, g2 behaves normally.

Does anyone have any idea about what happened?

Best,

Dan

[jrjoh...@gmail.com]

Hi

Getting a negative g2 is only possible when you get a significant occupation in the largest fock state in the truncated hilbert space, so yes I think you have used a too small space size. I would recommend plotting the fock state distribution for both modes for the steadystate density matrix. Then you could visually confirm whether your truncation is OK or not. If increasing the hilbert space size does not solve the problem it could be an issue with the steadystate solver.

Rob

[Paul Nation]

Another trick that I have used in the past to check if the number of states is enough is to look if the standard bosonic commutator for the two modes is close to one.  If you have a large occupation in the highest Fock state, then the commutator will be far from one indicating that truncation is an issue.

Regards,

Paul

[jrjoh...@gmail.com]

I just ran your program and had a look at the values of expect(ac.dag()*ac.dag()*ac*ac, rhoss) and expect(ac.dag()*ac, rhoss) and they are very close to zero both of them for the parameters in your attached example, so in this case I think the strange values of g2 is because of numerical errors (when dividing two very small numbers, ~1e-14 or so). 

rob

[dan hu]

Hi Robert and Paul, thank you very much!

Robert: For the numerical errors, is that possible it can change the sign of the dividing? Because I look at more details of the case of eta=0.01kappa and there are still cases where g2 is negative and the expect(ac.dag()*ac.dag()*ac*ac, rhoss) is also very small (shown 0 when I load it again, it might because when I save the data it truncate into 0).  Is there any thing else I can do to avoid this type of error besides increase the driving? Because if I increase the driving, I need to increase the subspace which will greatly slow down the code.

Paul: Do you mean to calculate \Delta = <a_c*a^\dagger_c> - <a^\dagger_c* a_c> and to check whether \Delta close to one?  

Best,

Dan

[dan hu]

Here are more details about the case where I got negative g2,  because expect(ac.dag()*ac.dag()*ac*ac, rhoss) = -3.3268308289763883e-13 and expect(ac.dag()*ac, rhoss) = 2.854220912679858e-07

the Fock distribution for the steady-state is in the attached figures,  where the fock distribution for the reduced cavity state has distribution for the negative number. What does it mean? Is it a sign of insufficient subspace for cavity field?

I checked the commutator relation for this steady-state and it is okay, expect(ac*ac.dag(), rhoss)-expect(ac.dag()*ac, rhoss) = 1.0000000000005131.

and expect(bm*bm.dag(), rhoss)-expect(bm.dag()*bm, rhoss) = 0.9967417881186798

From these information, can you tell me which subspace is insufficient?  

Best,

Dan

[Paul Nation]

The cavity does not have a negative number distribution.  You are using a bar plot and the width of your bars is too large.  Hence the ground state probability bar overflows into the negative.  The default width=0.8 which is exactly what you see.  Use “width=0.1” or so  in the call to bar() and you will get a better plot.

Since both are near unity, you are most likely not running into truncation issues, as your plots also indicate.  As Rob said, the issues you are dealing with are probably arising from numerical error.

Paul

--
You received this message because you are subscribed to the Google Groups "qutip" group.
To unsubscribe from this group and stop receiving emails from it, send an email to qutip+un...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
<CavityFock.png><MechanicalFock.png>

[dan hu]

Hi Paul,

Do you know for calculating g2, which order of  expect(ac.dag()*ac.dag()*ac*ac, rhoss ) and expect(ac.dag()*ac)**2 are reliable? 

When I increase the driving strength, the negativity of the g2 will disappear but I am not sure whether these results are still suffering from the numerical errors or not. Because in my code, I want to use as small driving as possible so that I could limit my subspace.

Best,

Dan

[jrjoh...@gmail.com]

Hi Dan

The problem seems to be that the limit of expect(ac.dag()*ac.dag()*ac*ac, rhoss)/expect(ac.dag()*ac, rhoss)**2 when expect(ac.dag()*ac, rhoss) goes to zero is not well defined numerically, even though it can be well defined analytically. You need to drive your system to at least have a small probability of being excited. As it is now in your calculation, the probability of being in the zero-photon state is:

 1.336599401472816e-09

so basically only the zero-photon fock state is occupied, but the operator ac.dag()*ac.dag()*ac*ac has a zero at this position, so when you do expect(ac.dag()*ac.dag()*ac*ac, rhoss) you end up summing a bunch of very small numbers. Because of numerical errors in the steadystate solver some of these can be negative. To see this you can do something like

rhoss.diag()

The diagonal of rhoss should of course be strictly positive, but as you can see there are some very small negative numbers due to numerical errors. This is usually not a problem but in your case with zero photon in the mode you are trying to calculate g2 for it gives you negative results. 

The only reasonable work-around that I can see is that you drive the system so that not only the zero-photon fock state is occupied (a very weak driving should be enough, so you should not have to increase the hilbert space because of this).

rob

[dan hu]

Hi Robert,

Thanks for the suggestions. I did several tests and it shows that  for the stronger driving, such as epsilon=0.1kappa, I have to increase the cavity subspace to nc=4. Because for any driving, the steady-state g2 must be close to 1 for zero optomechanical coupling g0/wm=0 ,  with the driving epsilon=0.1kappa, only when the cavity subspace is chosen as nc=4, the g2 = 0.998524. For nc=3, the corresponding g2=0.945532.  

For the subspace (nc, nm)=(4, 160), the steady-state with direct method run really slowly (about 1 hour ) and use a lot of memory ( around 6GB ). In this case, does the iterative method help?  Which iterative method do you suggest to use for calculation of steady-state g2?

Best,

Dan

[Paul Nation]

Your runtime is large because trying to find the LU factors of a matrix with 10^11 elements is computationally demanding.  Iterative methods are very tricky to use because the Liouvillian is non-symmetric and therefore all of the mathematical theorems involving iterative methods are not valid.  Therefore, I do not recommend this.  I am working on a different method that makes use of the Hermitian symmetry of the density matrix but this is not yet ready, and there is no guarantee that it would work any better.  In contrast to naive expectation, solving for the steady state solution is non-trivial.

In your case, you might be able to lower the 160 states of around 140 without too much loss in terms of accuracy.  This may help the runtime.  I have also found that umfpack tends to run faster than the SuperLU library used by SciPy, but then you would need to install scikits.umpfack yourself.

Paul

[Ka Wa Yip]

Hi Paul,

Sorry to bother you. I have encountered the error using the default steady-state solver (direct LU factorization). It is (RuntimeError: Factor is exactly singular). I changed the steady state method by adding( method='power', use_rcm=True). It does not have error but it produces strange graph. So I guess it would affect gmres or svd.

I just wonder how to deal with this error. Can you give me some hints on how to deal with the singular factor error in qutip by changing the hamiltonian or c_op_list? For example, adding noise to the singular matrix? 

Thank you very much for giving me advice. Thank you for the improvement in qutip3. 

[Paul Nation]

You should show us your code in order to figure out what is going on.  Your matrix should not be singular though.

[Myung-Joong Hwang]

Hi all, 

i ran into similar problems before, that is, the excited state occupation is too small that correlation functions becomes very sensitive to numerical errors. 

Can we get around this problem by using dtype with a higher precision than double? I guess it would be very costly, but I wonder if it is possible to choose custom dtype in the qutip and if it would help to get around these type of problems without changing physical conditions like driving strength.

myungjoong

[jrjoh...@gmail.com]

There is no easy way to change the dtype of the sparse matrices used in QuTiP. We already use np.complex128, so I don't think there is much to be gained even if you could easily change the dtype.​ I think the only solution is to include some driving so that at least you have a small nonzero occupation in the cavity, otherwise calculating the g2 correlation results in dividing two very small number, ~0.0 / ~0.0, which is not numerically robust.

Rob

[Paul Nation]

As Rob just mentioned, there is no easy way to do this.  Especially because when calculating the steady state of a system all methods in QuTiP call external libraries that have no way of handling multi-precision numbers.

Paul

-- 

----------------------------------------
Paul D. Nation
Assistant Professor
Korea University
Department of Physics
Anam-dong 5, Seongbuk-gu
Seoul 136-713, South Korea
Email: pna...@korea.ac.kr
Phone: +82-02-3290-3092
Web: http://nqdl.korea.ac.kr
----------------------------------------

[Akimasa Ihara]

Hi, I am trying to simulate the excitation spectrum of Ba-137 using qutip's steadystate function and am encountering the same issue where I get the error - RuntimeError: Factor is exactly singular. Adding method="power" works sometimes but gives weird results too. I've checked the determinant of the matrix and it is non-zero (though it is massive - on the order of 10^300ish since it is a 32x32 matrix). Did you guys ever figure out what the cause of the error was? Thank you so much for the help. 

[Simon Cross]

Hi Akimasa,

> Hi, I am trying to simulate the excitation spectrum of Ba-137 using qutip's steadystate function and am encountering the same issue where I get the error -
> RuntimeError: Factor is exactly singular. Adding method="power" works sometimes but gives weird results too.

First, I wanted to check that you are using a more recent version of
QuTiP than 3.0.0 -- i.e. the version mentioned in the subject of this
thread.

In general, as Paul mentioned, calculating these steady states is
numerically tricky, and that is one possible underlying cause of the
error you're seeing. It is also possible that there is a bug in QuTiP
or in the underlying solver being used, but it's impossible to say for
sure without a small piece of code that reproduces the problem and
some insight into what the correct output should be.

Yours sincerely,
Simon Cross