# Algorithm for Complete Eigenvalue Problem of a Real Symmetric nxn Matrix

I have a nxn covariance matrix (so, real, symmetric, dense, nxn). 'n' may be very very very big! I'd like to solve complete eigenvalue (+eigenvectors) problem for this matrix. Could somebody tell me what the fastest algorithm to do it?

P.S. I'd like to make GPGPU implementation using OpenCL. Typical sizes is 10000x10000 or even bigger. That's why I'm looking for fastest algorithm

P.S.S. I'd like to write own implementation (just for education) so please don't advise already existing libraries

• Could you give a typical size of 'n'? Is their a particular programming language that you want to use?
– brm
Jan 18 '14 at 11:37
• I'd like to make GPGPU implementation using OpenCL. Typical sizes is 10000x10000 or even bigger. That's why I'm looking for fastest algorithm. Jan 18 '14 at 11:45
• You might as well look into CULA if you are fine with a CUDA implementation. CULA : culatools.com Jan 18 '14 at 22:20

10000 by 10000 isn't really very big. A double precision matrix of this size requires 800 megabytes of memory. You'll need at least as much memory to hold the matrix of resulting eigenvectors, and you'll need additional working storage, but this is well within the capability of a typical desktop machine with 8 gigabytes or more of memory. It would also fit into the graphics memory on a higher end GPGPU.

In terms of time, this also isn't too bad. I just computed the eigenvalues and eigenvectors for a random symmetric 100000 by 10000 matrix on a desktop machine (using MATLAB/LAPACK and running with 8 cores), and it took about 3 minutes to run.

If you only need to do this computation for a small number of matrices of this size, then it's likely not worth your time and effort to move this to a GPGPU. However, if you plan on doing these computations 24/7, then putting some effort into speeding them up might be appropriate.

This is the sort of computation that you might be able to speed up with a high end GPGPU, but don't be surprised if it doesn't perform that well on a more basic graphics card. The issues here are that your graphics card might not have sufficient memory, and that double precision performance is not a strong point of many GPGPU's.

Since eigenvalue computation is commonly available through the LAPACK libraries, I'd suggest that you look for an implementation of the LAPACK libraries that has been optimized for your GPGPU. In particular, check out the OpenCL version of the MAGMA library:

Magma

• Thank you very much! But anyway... Could somebody tell me what the fastest algorithm for eigenproblem in my case(very big real symmetric matrix)? I'd like to write own implementation (just for education). Jan 20 '14 at 6:17

Regarding your question about the fastest algorithm for the complete symmetric eigenproblem...

There are two main phases to the computation. First is reduction to tridiagonal form using (blocked) Householder reflections. This runs in cubic time w.r.t the matrix dimension, and you can easily code up your own non-blocked version, while the blocked version is somewhat more complicated.

Computing the eigenvalues of the result tridiagonal matrix can be done in pseudo-quadratic time using the most cutting edge algorithms (MRRR or Divide and Conquer), or using the traditional QR or QL iteration. The tricky part is that you have to not only do the iterations (which is "trivial"), but you also have to deal with deflation and various other optimizations that can easily speed things up by an order of magnitude or more. I've only linked to (equivalent) Lapack code for the traditional stuff, the code for MRRR is incredibly subtle and can often fail, and the step in Divide and Conquer for solving the secular equation to glue the divisions together is also incredibly subtle. Coding these things up is highly nontrivial even for ordinary CPU execution; I can't imagine how hard it would be on a GPU.