Version 5.0.1 performance issue

[Katya Fouka]

Dear QuTiP community,

I have been curious to try the latest version of QuTiP but I am hesitant due to performance issues. My code seems to run much slower when using v.5.0.1 in comparison to v.4.7.5. For example, when studying the time evolution of laser cooling of a 3 level system coupled to 15 motional degrees of freedom, it requires 1 minute with v.4.7.5 and 25 minutes with v.5.0.1. Essentially the only change I make to the code is

#version 4.7.5

options = Options(nsteps=10000)
rho = mesolve(H,rho0,tlist,c_ops=c_ops,e_ops=[a.dag()*a],options=options,progress_bar=True)

#version 5.1.0

rho = mesolve(H,rho0,tlist,c_ops=c_ops,e_ops=[a.dag()*a],options={"progress_bar": "tqdm","nsteps": 10000})

I run the code on Windows using Anaconda and qutip.about() yields:

QuTiP Version:      4.7.5
Numpy Version:      1.26.4
Scipy Version:      1.11.4
Cython Version:     3.0.11
Matplotlib Version: 3.8.0
Python Version:     3.11.11
Number of CPUs:     12
BLAS Info:          Generic
OPENMP Installed:   False
INTEL MKL Ext:      False
Platform Info:      Windows (AMD64)

when using v.5.0.1 I keep the same versions for the other packages.

For reference, I am attaching my code here. Any leads will be much appreciated.

Thank you in advance,

Katya Fouka

[Katya Fouka]

Correction: where I write v.5.0.1 I mean v.5.1.0.

[Éric Giguère]

qutip 5.X can support many format for Qobj whereas v4 support only sparse arrays. It's not always great at using the right format.

What happen in your case is that it uses a Dense matrix for the Liouvillian where your matrix is quite sparse.

Unless expressly told, `Qobj.expm` will return a dense matrix. v4 will round off the small off-diagonal values automatically.

To get a sparse matrix, you can use (1j*y).expm(dtype="CSR") or (1j*y).expm().to("CSR").

To control the cutoff you can use tidyup: (1j*y).expm().to("CSR").tidyup(atol=1-e8), otherwise the cut off comes from floating point precision.

Another issue that v5 is bad with is creating operators from kets.

basis(3, 0) will use a Dense representation per default and this will spread to the whole Liouvillian.

Every Qobj creating function should have a dtype parameter to tell it the type you want:

e = tensor(basis(Natomic, 0, dtype="csr"), qeye(Nmotion))

r = tensor(basis(Natomic, 1, dtype="csr"), qeye(Nmotion))
g = tensor(basis(Natomic, 2, dtype="csr"), qeye(Nmotion))

In you case, you can convert everything once at the end:

H = H.to("CSR").tidyup(atol)

c_ops = [c2.to("CSR"), c3.to("CSR")]

To get timings that are within 10% of v4's.

Hope this helps

[Katya Fouka]

Hi Éric,

Thank you so much for your answer! Your suggestion seems to resolve the issue :)

Best,

Katya