How can I get eigenvalues and eigenvectors fast and accurate?

Question

I need to compute the eigenvalues and eigenvectors of a big matrix (about 1000*1000 or even more). Matlab works very fast but it does not guaranty accuracy. I need this to be pretty accurate (about 1e-06 error is ok) and within a reasonable time (an hour or two is ok).

My matrix is symmetric and pretty sparse. The exact values are: ones on the diagonal, and on the diagonal below the main diagonal, and on the diagonal above it. Example:

How can I do this? C++ is the most convenient to me.

Answer 1

Your system is tridiagonal and a (symmetric) Toeplitz matrix. I'd guess that eigen and Matlab's eig have special cases to handle such matrices. There is a closed-form solution for the eigenvalues in this case (reference (PDF)). In Matlab for your matrix this is simply:

n = size(A,1);
k = (1:n).';
v = 1-2*cos(pi*k./(n+1));

This can be further optimized by noting that the eigenvalues are centered about 1 and thus only half of them need to be computed:

n = size(A,1);
if mod(n,2) == 0
    k = (1:n/2).';
    u = 2*cos(pi*k./(n+1));
    v = 1+[u;-u];
else
    k = (1:(n-1)/2).';
    u = 2*cos(pi*k./(n+1));
    v = 1+[u;0;-u];
end

I'm not sure how you're going to get more fast and accurate than that (other than performing a refinement step using the eigenvectors and optimization) with simple code. The above should be able to translated to C++ very easily (or use Matlab's codgen to generate C/C++ code that uses this or eig). However, your matrix is still ill-conditioned. Just remember that estimates of accuracy are worst case.