I'm trying to find both the dominant $k$ left and right eigenvectors, that is,
$$V_L\mathcal{A} = \Lambda V_L\\ \mathcal{A}V_R = V_R\Lambda\\ V_LV_R = I_{k\times k}$$
- $V_L$ being the $k\times N$ matrix of left eigenvectors,
- $V_R$ being the $N\times k$ matrix of right eigenvectors, and
- $\Lambda$ being the diagonal $k\times k$ matrix of eigenvalues.
I'm using ARPACK (via scipy.sparse.linalg
). scipy.sparse.eigs
gives me the right eigensystem $V_R$ and $\Lambda$; I've tried two methods to get $V_L$:
-
- Apply
scipy.sparse.linalg.eigs
to $\mathcal{A}^T$ to get left eigenvectors and eigenvalues (call them $V_{L,0}$ and $\Lambda_L$) - If $V_{L,0}V_R$ is singular or $\Lambda_L$ is too different from $\Lambda$, reduce $k$; if that doesn't work, restart
- Set $V_L \leftarrow (V_{L,0}V_R)^{-1}V_{L,0}$; or
- Apply
- For each of the $k$ eigenvalues $\lambda_i$, apply
scipy.sparse.linalg.lsmr
to the overdetermined system $$\lambda_i^{-1}v_L^{(i)}\mathcal{A} = v_L,\\ v_L^{(i)}v_R^{(j)} = \delta_{i,j}$$ This is much more reliable than the other method, but it has a much worse best-case running time and it also has a quite limited accuracy.
I imagine it should in principle be possible to get left and right eigenvectors via IRAM; is this available in some way on ARPACK? (As mentioned I've been using scipy.sparse.linalg
, but I'm amenable to accessing ARPACK more directly as a last resort). Is there some way to repackage this approximate eigendecomposition as something that ARPACK can handle (sort of like the $AA^\dagger$ trick for SVD)? Or is there some other way to get $V_L$, after having gotten $V_R$ and $\Lambda$, that's reasonably efficient, reliable, and accurate?
On the matrices I'm trying to diagonalize:
- Complex, non-Hermitian (also not symmetric, unitary, normal, etc.)
- Represented via
scipy.sparse.linalg.LinearOperator
, not technically a sparse matrix - Equally efficient to perform left- and right- vector multiplication
- Magnitude of eigenvalues should theoretically decay fast (between polynomially and exponentially); initial eigenvalues tend to decay even faster
- For an example see arXiv:1510.00689; operator is defined in Fig. 2 and eigenvalue problem is described in Fig. A.4.b.