Complex Eigenvalues using eig (Matlab)

Question

I wanted to find and plot the eigenvalues of large matrices (around1000x1000). 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. Tridiagonal Toeplitz

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,'.')

Any help would be appreciated.

Kirill · Answer

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.09times 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$.

Federico Poloni · Answer

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.

Complex Eigenvalues using eig (Matlab)

2 Answers

Add your own answers!

Ask a Question