Numerical eigenvalue problem of a superoperator I have a linear super-operator $\mathcal{D}: \mathcal{M}_N(\mathbb{C}) \rightarrow \mathcal{M}_N(\mathbb{C})$ which acts on a square matrix $\rho$ and that returns a square matrix $\mathcal{D}(\rho)$. I want to compute the eigenvalues $\lambda_i$ and eigenvectors $\rho_i$ of $\mathcal{D}$ such that
$$
\mathcal{D}(\rho_i) = \lambda_i \rho_i
$$
and I am only interested in the first few eigenvectors of lowest eigenvalue amplitude (i.e. $\rho_0$ such that $|\lambda_0| < |\lambda_i|$). Is there a known package (in Python, Matlab, C/C++) that could solve this problem? If not, is there a known algorithm which I could use to write this eigensolver?
For instance, the super-operator could be a Lindblad operator:
$$
\mathcal{D}[F](\rho) = F\rho F^\dagger - \frac{1}{2} \{ F^\dagger F, \rho \}
$$
with $F \in \mathcal{M}_N(\mathbb{C})$.
 A: $\mathcal{M}_{N}(\mathbb{C})$ is a linear vector space with $\text{dim}(\mathcal{M}_{N}(\mathbb{C}))=N^2$. Since the superoperator acting on this space is also linear, the eigenvalue problem for this operator reduces to diagonalizing the $N^2\times N ^2$ matrix representantion of $D$ in terms of a basis of the aforementioned space.
To express $D$ as an explicit matrix, first choose a basis of $\mathcal{M}_N(\mathbb{C})$,$\{\mathbf{e}_{i}, ~~~i=1,...,N^2\}$ and compute the following coefficients:
$$D[F]({\mathbf{e}_i})=\sum_j C_{ij}\mathbf{e}_j$$ 
and represent the superoperator as the $N^2\times N^2$ matrix
$$(\mathcal{D}_F)_{ij}=C_{ij}$$
For this matrix, the solution of the nonsingular eigenvalue problem has $N^2$ solutions
$$D_F \rho=\lambda \rho\iff(\lambda_i,\rho_i=\sum_{j=1}^{N^2}S_{ij}\mathbf{e}_j)~~~,~~~ i=1,...,N^2$$
with the $\rho_i$'s being the matrices solving the eigenvalue problem of the superoperator expressed in the basis chosen. 
This algorithm doesn't need to be reimplemented, all you need is to recycle the classic eigendecomposition algorithms that already exist and if you are interested in the smallest matrix elements, you can use the Arnoldi-Lanczos algorithm as usual.
For a concrete choice of basis, you can choose
$$(\mathbf{e}_{i(j,k)})_{j'k'}=\delta_{jj'}\delta_{kk'}~~~, ~~~j,j',k,k'=1,...,N~~~,~~~i(j,k)=N(k-1)+j $$
where $i(j,k)$ is a typical enumeration choice of the eigenvectors, and it can be easily shown to take all integer values in the interval $1\leq i\leq N^2$, as expected. For example, we can compute the matrix elements of the Lindblad operator in this basis for the result:
$$(\mathcal{D}_F\mathbf{e}_{i(j,k)})_{lm}=F^{\dagger}_{lj}F_{km}-\frac{1}{2}\Big[(F^{\dagger}F)_{lj}\delta_{mk}+\delta_{lj}(F^{\dagger}F)_{km}\Big]$$
which can be directly fed into the diagonalization algorithms. I hope this is helpful!
