Eigensystem returns vectors which are not eigenvectors

Question

Short synopsis: for a specific family of sparse matrices, the eigensolver seems to be unstable (kernel quitting) for certain examples, and when it works it seems to consistently return vectors which are not eigenvectors.

I would like to know if there are any known options or pre-conditioning methods which fix this, and which are feasible for large sparse matrices.

More details:

I have a family of sparse Hermitian cyclic-banded matrices $M$, and I want to calculate the smallest (absolute value) eigenvalue for each of them.

However, the kernel seems to unexpectedly quit (a problem which I have not diagnosed) during execution of Eigensystem.

In a possibly related problem, when Eigensystem executes without quitting the Kernel, it consistently returns vectors which are not eigenvectors.

An example $1220\times1220$ matrix which can be obtained using Import (not Get) can be found here. The matrix is hermitian (see In[4] in image below) converting to a dense matrix and solving yields an eigenvector $\phi$ (In[5]) which satisfies the eigenvector equation (Out[7]).

Working with sparse matrices is much faster (compare $\tau$ from In[5] and In[8]) but is not even close to satisfying the eigenvector equation (Out[10]).

Using Method->banded resolves the issue (Out[13]), but the performance is horrible ($\tau$ from In[11] is worse than the dense solve).

I believe the default method is Arnoldi, and I have been unable to fix this by playing around with the options. Arnoldi always returns a vector which does not satisfy the eigenvector equation, even if it is handed the correct solution as the sarting vector.

Is there anything I can do about this problem?

Answer 1

Perhaps you can add a "Shift" option to help out the Arnoldi algorithm:

{{λ}, {v}} = Eigensystem[
    M,
    -1,
    Method->{"Arnoldi", "Tolerance"->10^-12, "Shift"->0}
]; //AbsoluteTiming

{0.007198, Null}

I added the tolerance to get a slightly better result than the default. Check:

λ
Norm[M.v - λ v]

0.0039334

3.96868*10^-15