Complex Eigenvalues using eig (Matlab)

Question

I wanted to find and plot the eigenvalues of large (around $1000\times1000$) matrices. But discovered when using the eig function in matlab, it gives complex eigenvalues when it shouldn't. For example, in the code below I have a Tridiagonal Toeplitz matrix, which should have all real eigenvalues.

But it seems eig is unstable for n=90 and returns a small complex error in a few of the eigenvalues. Is there a way I can get the eigenvalues more accurately?

clear parameters
close all
clc
n=90;
dd=-2.*ones(n,1);
ud=1.8*ones(n,1);
ld=.1*ones(n,1);
A = spdiags([ld dd ud],-1:1,n,n);
C=full(A);
g=eig(C);
g=sort(g);
cond(C)
plot(g,'.')

Answer 1

You have an ill-conditioned eigenvalue problem. Consider a perturbation $\delta A$ in your matrix $A$—with double-precision floats this is $O(10^{-16})$. It turns the eigendecomposition from $$ X^{-1}A X = \Lambda $$ into (approximately keeping $X$ fixed) $$ X^{-1}(A+\delta A)X = \Lambda + \delta \Lambda. $$ This means that the change in eigenvalues is bounded by $$ \|\delta \Lambda\| \leq \|X^{-1}\| \|X\| \|\delta A\| = \kappa(X)\|\delta A\|, $$

where $\kappa$ is the condition number.

Assuming that the perturbation in $A$ is on the order of machine epsilon, which is about $10^{-16}$, there will be an error on the order of $\kappa(X) \times 10^{-16}$ in the eigenvalues. This is not instability in the numerical method of eig, but rather a property of your mathematical problem. For example, you wouldn't have this problem if $A$ was symmetric, in which case $\kappa(X)=1$.

In your case, I calculated the condition number $\kappa_2(X)$ with extra precision for $n=64$ ($n=90$ failed to converge), and already it is $4.09\times 10^{39}$, which is much too large for ordinary double-precision floating-point arithmetic.

If you really need to solve such ill-conditioned eigenvalue problems, the easiest way is probably to use enough extra precision to keep $\kappa(X)\epsilon_{\mathrm{mach}}$ small, so you'd need about 60 digits for $n=64$ and more for larger $n$.

Answer 2

If the superdiagonal and the subdiagonal have the same sign, like in your example, you can scale your matrix to be symmetric. Replace your matrix $C$ with $DCD^{-1}$, with $D=\operatorname{diag}(1,d,d^2,\dots,d^{n-1})$, choosing $d$ so that it becomes symmetric. (I believe it's going to transform $\operatorname{tridiag}(a,b,c)$ into $\operatorname{tridiag}(\sqrt{c/a},b,\sqrt{c/a})$, but I am too lazy to check.)

Then Matlab will use algorithms for the symmetric eigenproblem, which are faster, more accurate, and always return real eigenvalues.