Why are my eigenvalues complex? (eig)

I wanted to find and plot the eigenvalues of large matrices (around1000x1000). But discovered when using the eig function, it gives complex eigenvalues when it shouldn't. 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.

Answer 1

[Christine Tobler]

The eigenvalues of a real matrix are only real if the matrix is symmetric. The matrix C is not symmetric, therefore the eigenvalues are either real or complex conjugate pairs.

[I Ahmed]

But a tridiagonal Toeplitz matrix, according to the wikipedia page, should have real eigenvalues, if the upper and lower diagonals are the same sign. Tridiagonal-Toeplitz-Matrix

[Geert Awater]

That's correct, and for n up to 75 your code produces in fact real eigenvalues. Beyond that numerical noise results in complex eigenvalues.

To get the correct, real eigenvalues for (much) larger n, apply a similarity transform to A to make it symmetric whie keeping it tridiagonal (see the "tridiagonal matrix" lemma in wikipedia) and put the resulting sparse matrix directly into the eig function.

    n=10000;
    dd=-2.*ones(n,1);
    ud=1.8*ones(n,1);
    ld=.1*ones(n,1);
    A = spdiags([sqrt(ld.*ud) dd sqrt(ld.*ud)],-1:1,n,n);
    g=eig(A);
    condest(A)
    plot(g,'.')

[Bruno Luong]

Geert Awater comment deserves to be on a separate answer (and accepted)

[Bruno Luong]

The condition number of A is not relevant in eigenvalue computation, what is more relevant is the condition number of the eigen-vectors matrix.

When they are large; the eigen spaces are almost parallel and it causes numerical algoritm to fail. The similarity transformation is the excellent trick to improve the conditioning without changing the eigen-values.

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);

[V,g]=eig(full(A));
cond(V) % 1.4975e+53

D=spdiags([sqrt(ld.*ud) dd sqrt(ld.*ud)],-1:1,n,n);
[W, g]=eig(full(D));
cond(W) % 1.0000