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?

$\endgroup$
8
  • $\begingroup$ Not sure if relevant, but these matrices are obtained from blocks of a covariance matrix like this $\endgroup$ Commented Nov 13, 2020 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$ Commented Nov 14, 2020 at 0:43
  • $\begingroup$ Indeed, scipy.linalg.schur seems to be similar speed to scipy.linalg.eig $\endgroup$ Commented Nov 14, 2020 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$ Commented Nov 14, 2020 at 3:29
  • $\begingroup$ The scipy routine has an option "output='real'" to do the real Schur decomposition. $\endgroup$ Commented Nov 14, 2020 at 3:43

1 Answer 1

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.

$\endgroup$
2
  • $\begingroup$ thanks, nice trick! Cholesky was 10x faster than eigh and inverse 5x, so this reduction pays off in computation time $\endgroup$ Commented Nov 17, 2020 at 19:41
  • $\begingroup$ @YaroslavBulatov You make a call to chol and one to eigh with this solution. $\endgroup$ Commented Nov 17, 2020 at 19:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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