I have been trying to compute eigenvalues of a rather sizable matrix A
, about $500 \times 500$ (but sparse). I asked Mathematica to compute Eigenvalues[A]
, and left it to work. After a night of computation, Mathematica still failed to produce the answer. Just out of curiousity, I tried to see what happens if I replace A
by it's numeric approximation (before, A
had integer entries, hence exact). (I didn't quite expect it to make a difference, but I was trying some random things). Much to my surprise, the answer appeared immediately. Hence my question: How does it happen that Mathematica produces the answer so much faster for inexact input? To what extent can this output be trusted?
1 Answer
Calculating eigenvalues involves solving for the roots of the characteristic polynomial, which is of degree equal to the order of the size of the matrix. When you input real numbers, it can search for the roots of the polynomial using numerical techniques. When you input exact integers (or rationals, probably) it tries to find exact answers for the roots of the polynomials. This is a much more difficult task, since closed form solutions only exist up to order 5 or so.
In terms of trusting the answers, J.M.'s suggestions are good. It would also make sense to check the answers using alternative methods. For instance, the product of the eigenvalues is equal to the determinant. So you could take the answer you get for the eigenvalues and multiply them together; if they are close to the Det
then this provides some confidence. If they are not the same, then the answer can be discarded.
-
2$\begingroup$ Actually, for quintics and higher degrees, Mathematica is certainly forced to express the results in terms of
Root[]
, whereupon the effort is in determining the (coefficients of the) minimal polynomial corresponding to the eigenvalue, plus root isolation. All this takes a nontrivial amount of work. $\endgroup$ Commented May 7, 2013 at 11:45 -
2$\begingroup$ At least for inexact arithmetic, I'm pretty sure Mathematica never looks for the roots of a polynomial (because that's not numerically stable). It probably uses an iterative method like Jacobi's. $\endgroup$ Commented May 7, 2013 at 12:34
-
$\begingroup$ I remember taking a class on numerical methods. We read a paper called "18 Ways not to Calculate the Matrix Exponential" -- the various complexities of eigenvector calculations also played a large role. So yes, I'm sure they are doing all sorts of clever things. $\endgroup$– bill sCommented May 7, 2013 at 12:44
-
$\begingroup$ @nikie, at least for eigenvalues, it uses the QR algorithm in the inexact case. $\endgroup$ Commented May 7, 2013 at 12:46
-
$\begingroup$ Nice summary of the QR algorithm at wikipedia: en.wikipedia.org/wiki/QR_algorithm $\endgroup$– bill sCommented May 7, 2013 at 13:35
Eigenvalues[N[A, 20]]
), 30 digits, ... and so on. If there is some agreement between these results, then I trust the results a little bit more. $\endgroup$RootApproximant
may save you though. $\endgroup$