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
A = spdiags([ld dd ud],-1:1,n,n);
  • 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, 2017 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$ Sep 2, 2017 at 8:13

2 Answers 2


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

  • $\begingroup$ +1 that's a very valuable answer! I do have a follow-up question though: You say that the condition number of the eigenvector matrix is 1 if the matrix is symmetric, but eigenvectors can be multiplied by a constant and still be eigenvectors: A = [ 5 1i ; -1i 8]; [V D] = eig(A); cond(V); gives 1 but after V(:,1) = 10*V(:,1) it becomes 10. I suppose MATLAB by default normalizes the eigenvectors in V (your X) though. Also in your last equation, should it be $2\kappa$ since we have both $||X||$ and $||X^{-1}||$? $\endgroup$ Apr 13, 2021 at 22:12
  • $\begingroup$ Likewise cond(inv(V))=10 so it seems, so it's not like the factor of 10 cancels out with a factor of 1/10 in $X^{-1}$. I guess it would help if you told us which matrix norm $||\cdot|| $ you're using. $\endgroup$ Apr 13, 2021 at 22:34
  • 1
    $\begingroup$ @user1271772 : The eigenvalue algorithms usually are (extremely) evolved variants of the QR algorithm. All transformations inside the algorithm are orthogonal. As a symmetric matrix has an eigen-decomposition with orthogonal matrices, the algorithm will return this orthogonal transformation matrix as $X$. For non-symmetric matrices the decomposition with the orthogonal transformation results in an upper triangular matrix, $Q^TAQ=R$ (with 2x2 diagonal blocks for complex eigenvalues). Diagonalization of $R$ results in a non-orthogonal $X$. $\endgroup$ Apr 14, 2021 at 7:08
  • $\begingroup$ @user1271772 : (cont.) That root clusters around zero become near indistinguishable under floating-point noise from a near multiple root and thus mutate to near-regular stars in the complex plane is a property of the matrix $A$ and its condition. $\endgroup$ Apr 14, 2021 at 7:11

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.

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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