6

I am looking for an C/C++/Python algorithm implementation that calculates eigenvalues and eigenvectors of a symmetric, positive semidefinite covariance matrix.

A general-purpose eigen-decomposition algorithm has about $O(n^3)$ complexity, but maybe a faster method exists for symmetric, positive semidefinite covariance matrices.

Rodrigo de Azevedo
  • 1
  • 4
  • 13
aleksv
  • 201
  • 2
  • 7
  • 3
    A real matrix is a covariance matrix iff it is symmetric positive semidefinite. So you are asking for eigen-decomposition of a symmetric positive semidefinite matrix. – Mark L. Stone May 10 '18 at 20:54
  • Yes, that's seems to be the case. – aleksv May 10 '18 at 21:09
  • What's the size of your matrix? Is the matrix sparse? Do you need all eigenpairs or just some (which ones?). For dense symmetric matrices, LAPACK has efficient implementations of several state of the art algorithms, like divide and conquer, QR and bisection + inverse iteration. Some of those routines are specifically designed for SPSD matrices. See here. – GoHokies May 10 '18 at 21:35
  • I need all eigenpairs, matrix is not sparse. – aleksv May 10 '18 at 21:54
  • 2
    in my algorithm the covariance matrix is changed every iteration, presumably by a small value - if your updates are of small rank $k \ll N$, then the eigendecomposition can be updated with ${\cal O}(k N^2)$ FLOPs, see this q&a. For small-norm full-rank updates and $N \approx 100$ the fastest approach will likely be to re-calculate the eigendecomposition. – GoHokies May 11 '18 at 11:33
  • The updates are small norm full-rank. I'll probably simply reduce period of eigen-decomposition and use available O(N^3) methods. – aleksv May 11 '18 at 12:30

4 Answers4

3

TL;DR: I don't think you will be able to beat any real/vendor eigensolver, but it's fun to think about.

Deeper cut: One classic algorithm for symmetric eigendecomposition is tridiagonalizing A=QTQ' via householder methods, followed up by QR iteration upon T. The QR iteration is (very) loosely based on the iteration: [Q,R] = A; A = R*Q. That is, alternating between QR decomposition and then multiplying them out in reverse order. In the limit of many iterations, A will converge to a diagonal matrix (thus displaying the eigenvalues) and is also similar (same eigenvalues) to the original input.

For symmetric positive definite A, I think you could in theory beat this algorithm using a treppeniteration-like method based on Cholesky decomposition [Consult Golub & Van Loan 3rd ed, chapter 8, problem 8.2.1]. It would be (very) loosely based on the iteration [G,G'] = chol(A); A = G'*G. That is, computing a Cholesky decomposition and then multiplying them in reverse order. Remarkably, this converges to a diagonal matrix too, which is similar to the original input. The departure from orthogonal iterations is mild cause for concern, but fortunately the Cholesky decomposition is quite stable, too. You would also want to "frontend" this algorithm using householder tridiagonalization, so that all the A's in question become tridiagonal and the Cholesky's are all banded ones, with band=1. Fundamentally, band=1 Cholesky is less flops/simpler than band=1 QR via givens rotations, so that's how you could (possibly) come out on top.

I think you could also "accumulate" all these similarity transforms as you go (essentially, banded backsolution steps), to build the eigenvectors. This is much how the givens rotations in are accumulated in the classic QR iteration. If this is unworkable for some reason, you can always use inverse iteration at the end once you have the eigenvalues in hand.

All that said, real/vendor solvers are not just performing the dumb QR/RQ/QR/RQ iteration .. there's (at the very least) shifting, plus a whole wider class of (faster!) algorithms in the field (divide and conquer, MRR, bisection, etc). This is a ridiculous amount of machinery/cumulative improvements to try to compete against .. many (!!) man-years of effort have been put into the development of these algorithms and their implementations (EISPACK/LAPACK/MKL/etc). I think it's an interesting thought experiment (wow, an eigensolver built from such a simple decomposition) but not a very practical one.

rchilton1980
  • 4,862
  • 13
  • 22
  • Thanks for input. I'm using routines tred2 (Householder tridiagonalization) and tql2 (QL algorithm) from Jama package. But they run at O(N^3) complexity which is a bit steep for my needs. Your Cholesky decomposition suggestion seems interesting, but as it requires tridiagonalization, it will still be around O(N^3). The problem, my covariance matrix may be 100x100 and I need to evaluate a lot of such matrices. – aleksv May 10 '18 at 22:43
  • 1
    @aleksv How many 100 by 100 matrices? I can bang out almost 1000 full eigendecompositoions/sec using MATLAB eig on a vanilla PC. – Mark L. Stone May 10 '18 at 22:58
  • Probably, tens of thousands in the run of the "outer" algorithm that makes use of eigenpairs. – aleksv May 10 '18 at 23:01
  • @alexv, yep, all these algorithms are O(N^3) – rchilton1980 May 10 '18 at 23:28
  • @aleksv at $100\times 100$ size, you probably should not worry about $O(N^3)$, but about the constants and lower order terms that might dominate at such small sizes. – Anton Menshov May 11 '18 at 02:30
2

There are specialized methods for the eigendecomposition of symmetric matrices. LAPACK has DSYEV, Numpy has numpy.linalg.eigh. They are still $O(n^3)$, but overall cheaper and more accurate. You should use them.

I am not familiar with JAMA, but from what I could Google and understand in a minute it is a pure-java port of an old pre-LAPACK implementation of the QR/Francis algorithm. That might be what is slowing you down. Try replacing it with native code: find whatever library wraps LAPACK in your language, and use it.

As far as I know, you can't get any further speedup from the fact that your matrix is positive semidefinite nor from the fact that it has diagonal 1 (a property that is sometimes assumed when people say "covariance matrix"). If you could, cheap scaling and shifting tricks would permit to apply the same speedups to all symmetric matrices, I think.

Federico Poloni
  • 11,344
  • 1
  • 31
  • 59
0

I think, for $100 \times 100$, you should be happy with the eig method in Matlab, or numpy.linalg.eig in Numpy/Python. Since it is a covariance operator, it is symmetric positive semidefinite. Probably you can speed up your implementation by passing this information through the solver.

For such small matrices there is little chance to solve it faster.

There are indeed methods to solve the problem faster, if your covariance matrix has some structure - and is much bigger than $100 \times 100$:

  • If your matrix is Toeplitz or Block-Toeplitz with Toeplitz blocks, circulant embedding can be used to sample in $O(N \log N)$. Link to a tutorial.

  • There are multipole methods that can be used to speed-up eigendecompositions of covariance operators. See this technical report.

  • You need to solve the eigendecomposition for several matrices. If there is some underlying structure, you could solve the problem on a reduced basis. There is a preprint from me and my colleagues that may help you.

Jonas
  • 134
  • 5
  • Interesting. Actually, in my algorithm the covariance matrix is changed every iteration, presumably by a small value, so on every iteration I have a matrix of differences between previous and new covariance matrix state. Maybe a method is available that may efficiently recalculate the eigenpairs based on the covariance matrix deltas. – aleksv May 11 '18 at 09:12
-2

The Schur decomposition is a good solution for finding the eigen values and eigen vectors of a symmetric matrix. In MATLAB, use schur.m.

AJ Johnson
  • 1