2
$\begingroup$

I have non-symmetric real-valued matrices with real-valued eigensystems. How to compute eigenvectors efficiently?

Using scipy.linalg.eig (which calls dgeev) is 3-4 times slower than scipy.linalg.eigh (=>dsyevd) and 2x slower than scipy.linalg.svd (=>gesdd), but neither SVD nor hermitian eigenvalue decomposition are appropriate here. I'm guessing the need to support complex valued storage is responsible for extra overhead. Is there a faster routine to use for matrices which have real-valued eigenvalue decomposition?

Improve this question
$\endgroup$
  • $\begingroup$ Not sure if relevant, but these matrices are obtained from blocks of a covariance matrix like this $\endgroup$ – Yaroslav Bulatov Nov 13 '20 at 22:29
  • $\begingroup$ You might try finding the real Schur decomposition of your matrix and looking to the diagonal for the eigenvalues, but I'm afraid that you'll find apparent complex eigenvalues (showing up as 2x2 blocks) due to round-off error. This would avoid storing complex matrices, but might not be any faster. $\endgroup$ – Brian Borchers Nov 14 '20 at 0:43
  • $\begingroup$ Indeed, scipy.linalg.schur seems to be similar speed to scipy.linalg.eig $\endgroup$ – Yaroslav Bulatov Nov 14 '20 at 0:58
  • $\begingroup$ Note that Schur decomposition routines can produce a complex or real result (In MATLAB, schur(A,'real') vs. schur(A,'complex') Using the 'real" option might make it a bit faster. $\endgroup$ – Brian Borchers Nov 14 '20 at 3:29
  • $\begingroup$ The scipy routine has an option "output='real'" to do the real Schur decomposition. $\endgroup$ – Brian Borchers Nov 14 '20 at 3:43
6
$\begingroup$

The trick is trying to find out why that matrix has real eigenvalues in the first place. Usually it is because a suitable set of conjugations turns it into a symmetric matrix, and then you can reduce to a symmetric computation.

Multiplying and dividing by $(AC)^{-1}$ you can rewrite $$ D_1 = C^{-1}BAC (C^{-1}BAC+I)^{-1}, $$ so your computation is equivalent to finding the eigenvalues $\mu_i$ of $M = C^{-1}BAC$ and then computing $\lambda_i = \frac{\mu_i}{\mu_i + 1}$ for each of them. Clearly the eigenvalues of $M$ are equal to those of $BA$. If $A=LL^T$ is a Cholesky factorization then those are equal to the eigenvalues of $L^TBL$, which is symmetric. Bingo.

For eigenvectors, the same computation shows that if $v$ is an eigenvector of $L^TBL$ then $(L^TC)^{-1}v$ is an eigenvector of $M$, and hence of $D_1$.

EDIT: oops, one does need an inverse in the last part.

Improve this answer
$\endgroup$
  • $\begingroup$ thanks, nice trick! Cholesky was 10x faster than eigh and inverse 5x, so this reduction pays off in computation time $\endgroup$ – Yaroslav Bulatov Nov 17 '20 at 19:41
  • $\begingroup$ @YaroslavBulatov You make a call to chol and one to eigh with this solution. $\endgroup$ – Federico Poloni Nov 17 '20 at 19:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.