I have a complex, non-Hermitian matrix $\mathbf{A}$, for which I need to find a few eigenvalues and eigenvectors in the generalised eigenvalue problem:
$$\mathbf{A}\cdot \mathbf{x} = \lambda \mathbf{B} \cdot \mathbf{x},$$
with real symmetric positive definite matrix $\mathbf{B}$. I am currently using scipy.sparse.linalg.eigs, which is a Python wrapper around ARPACK routines.
I have successfully found both the eigenvalues and right eigenvectors, however I also need the corresponding left eigenvectors satisfying:
$$ \mathbf{y}\cdot\mathbf{A} = \lambda \mathbf{y}\cdot \mathbf{B}.$$
(According to my notation there is no complex-conjugation required on these left eigenvectors, but obviously undoing a complex-conjugate is easy enough).
As far as I can see, neither scipy nor ARPACK provides a ready-made solution to this problem, not could they be easily modified to give both sets of vectors. The possible solutions I have come across are:
Transpose my matrix $\mathbf{A}$ and make another call to
eigs. This is obviously inefficient, and I'm not sure if the returned eigenvalues are guaranteed to have the same values and be in the same order, so matching the left and right eigenvectors may not be reliable.Use inverse iteration, extrapolating from Numerical Recipes 2nd ed Section 11.7. Solve
$$(\mathbf{A}^T - \lambda\mathbf{B})\cdot \mathbf{y} = \mathbf{b}$$
iteratively for $\mathbf{x}$, starting from some random vector $\mathbf{b}$. The problem with this approach is that it may not find a full set of eigenvectors in the case of degenerate eigenvalues, which I expect to see in some cases.
Are there any better solutions available, or ways that my proposed solutions could be fixed to overcome their (perceived) limitations?
