scipy's sparse eigs yielding wrong eigenvectors

Question

I am a bit confused regarding the following issue:

I am computing fixed points of a quantum channel, which means I want to compute the leading eigenvector of a specific matrix. The matrix is such that its dimensionality is n^2 x n^2 and defined in such a way that the leading eigenvalue reshaped to a matrix with shape n x n is a positive matrix (self adjoint with positive eigenvalues).

If I do this with scipy.sparse.linalg.eigs however I get wrong results. The exact computation (using scipy.linalg.eig) however works fine. I tried around playing with the arguments for k and ncv for the solver, but didn't get i working properly unless I set k=n**2 in which case eigs just refers to eig. This, however, won't work in the case that I have actually in mind where the channel (super_op in the script below) is actually encoded as a LinearOperator. So I rely on using eigs :/

Anybody any idea how to get this run properly?

Thanks to everybody in advance!

import numpy as np
from numpy.random import rand
from numpy import tensordot as td
from scipy.sparse.linalg import eigs
from scipy.linalg import eig


n = 16
d = 3

kraus_op = .5 - rand(n, d, n) + 1j * (.5 - rand(n, d, n))
super_op = td(kraus_op, kraus_op.conj(), [[1], [1]]).transpose(0, 2, 1, 3)

########
# Sparse
########

vals, vecs = eigs(super_op.reshape(n**2, n**2), k=n*(n-1), which='LM')

rho = vecs[:,0].reshape(n, n)

print('is self adjoint: ', np.allclose(rho, rho.conj().T))

super_op_times_rho = td(super_op, rho, [[2, 3], [0, 1]])

print('super_op(rho) == lambda * rho :', np.allclose(rho, super_op_times_rho/vals[0]))

########
# Exact
########

vals, vecs = eig(super_op.reshape(n**2, n**2))

rho = vecs[:,0].reshape(n, n)

print('is self adjoint: ', np.allclose(rho, rho.conj().T))

super_op_times_rho = td(super_op, rho, [[2, 3], [0, 1]])

print('super_op(rho) == lambda * rho :', np.allclose(rho, super_op_times_rho/vals[0]))

the result is:

is self adjoint:  False
super_op(rho) == lambda * rho : True
is self adjoint:  True
super_op(rho) == lambda * rho : True

For completeness:

Python 3.5.2 numpy 1.16.1
scipy 1.2.1

Answer 1

After all I found the solution with some help of my colleagues:

While eig gives the eigenvectors (1) sorted (by magnitude) and (2) such that the first entry is real, eigs has another sorting that must not be as in eig and also does not regularize the complex phase of the eigenvector. Correcting the phase is easily done by dividing the tensor by the phase of the first entry (corresponding to ensuring that the first diagonal element is real to get rid of the freedom of choosing the complex phase of an eigenvector and making Hermiticity possible...).

So the corrected code snippet for the sparse case is:

vals, vecs = eigs(super_op.reshape(chi**2, chi**2), k=chi*(chi-1), which='LM')
# find the index corresponding to the largest eigenvalue
arg = np.argmax(np.abs(vals))
rho = vecs[:,arg].reshape(chi, chi)
# regularize the output array 
rho *= np.abs(rho[0, 0])/rho[0, 0]