3
$\begingroup$

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.

$\endgroup$
  • 3
    $\begingroup$ If you really are studying a tridiagonal Toeplitz matrix, you probably already know that its characteristic polynomial is expressible in terms of Chebyshev polynomials, and thus the roots are expressible as trigonometric functions of angles. Right? $\endgroup$ – J. M. Sep 1 '17 at 23:23
  • 3
    $\begingroup$ By the way, you know that cond(C) has nothing to do with the conditioning/sensitivity/stability of the eigenvalue problem, right? $\endgroup$ – Federico Poloni Sep 2 '17 at 8:13
6
$\begingroup$

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\|. $$ 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$.

$\endgroup$
4
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ This is what the routine FIGI() did in EISPACK. $\endgroup$ – J. M. Sep 2 '17 at 11:24

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.