mabuchilab / QNET

Computer algebra package for quantum mechanics and photonic quantum networks
https://qnet.readthedocs.io/
MIT License
71 stars 23 forks source link

SLH_to_qutip discards Lindblad operators #57

Closed ngedwin98 closed 7 years ago

ngedwin98 commented 7 years ago

Certain nontrivial Lindblad operators are discarded by SLH_to_qutip. An example is

a = Destroy(hs=LocalSpace("fock",dimension=10))
slh = SLH([[1]], [sympy.sqrt(2)*a], 0)
H, Ls = SLH_to_qutip(slh)

The Hamiltonian seems fine, but Ls is the empty list [], which is incorrect.

The bug appears to be caused by lines 200--201 in convert/to_qutip:

if L_qutip.norm() > 0:
    Ls.append(L_qutip)

QuTiP evaluates the trace norm for operators, and since sp.sqrt(2)*a has zero trace, it is not included in the returned numerical Ls.

goerz commented 7 years ago

Wait, what? Shouldn't norm(a) = 0 if and only if a=0 for any norm, whether Frobenius (trace) or one-norm, or anything else? Isn't the trace norm of L equal to sqrt(tr(L^\dagger L))? Is this a bug in qutip?

ngedwin98 commented 7 years ago

Hmm, yes you're right. I take my previous comment back -- the trace norm should be (L.dag()*L).sqrtm().tr(), and the fact that L has zero trace is not relevant.

However, at least with the version of QuTiP that I'm using (4.1.0), qutip.destroy(10).norm() produces 0.0, which is incorrect. This could be a bug in QuTiP -- I'll look into it a bit more.

goerz commented 7 years ago

We can switch to another norm to work around this, of course

ngedwin98 commented 7 years ago

So it seems like QuTiP is indeed using the wrong formula for the trace norm. In lines 956--958 of qobj.py, the trace norm is defined as

vals = sp_eigs(self.data, self.isherm, vecs=False,
    sparse=sparse, tol=tol, maxiter=maxiter)
return np.sum(sqrt(abs(vals) ** 2))

However, for a general non-Hermitian matrix, the trace norm is the sum of the singular values, not the magnitudes of the eigenvalues. This formula only generally holds for Hermitian matrices. To find the singular values of L, the correct procedure is to find the square roots of the eigenvalues of L.dag()*L, which is the formula we discussed above.

I can look into raising this issue on QuTiP's side. In the meantime, should QNET use a different norm or some other workaround, or wait for a fix from QuTiP?

goerz commented 7 years ago

The fix for this is scheduled to be part of the v4.3 release of qutip, see https://github.com/qutip/qutip/issues/748.

Nonetheless, I will switch the QNET code to use the 'max' norm, for compatibility with the current (and older) versions of qutip.