Let $M$ be a matrix which has the following properties:

1) $M$ is Hermitian 2) $M$ has only rational entries 3) $M$ is known to have rational eigenvalues

What algorithms are there for exactly computing the eigenvalues of $M$?

Mathematica seems to be able to diagonalise such matrices very very quickly, so such methods must exist (I am hoping they are not trade secrets of Wolfram Inc!). I am sceptical that Mathematica is using the rational roots test on the characteristic equation since the numbers you'd need to factorise get large quickly.

  • $\begingroup$ How big a matrix does it work for (in Mathematica)? $\endgroup$
    – Victor Liu
    Sep 12, 2013 at 18:20

2 Answers 2


The characteristic polynomial of a matrix $M$ can be written as

$\chi_M(z) = z^n + \textrm{trace(M)}z^{n-1}+\ldots + \det(M)$.

Since you know ahead of time that all the entries of $M$ are rational, you can be assured that $\det(M)$ is rational; say $\det(M) = p/q$. Now, the rational roots theorem assures you that, provided $\chi_M$ has a rational root $\lambda = k/l$, then $k$ divides $p$ and $l$ divides $q$. But, you know ahead of time (somehow) that all eigenvalues of $M$ are rational.

You can compute the determinant of $M$ much faster than the full characteristic polynomial by computing the LU-factorization of $M$. Since all the entries are rational, the entries of the triangular factors are also rational; you would of course have to ensure that the factorization is done symbolically so that the determinant is also computed exactly.

Trying all rational numbers $k/l$ such that $k | p$ and $l | q$ is also prohibitively expensive. However, you could use a conventional eigensolver to get a set of approximate eigenvalues $\hat{\lambda}_1,\ldots,\hat{\lambda}_n$, then look for the exact, rational eigenvalues with the desired divisibility properties nearby to each approximate one.

A cursory Google search didn't turn up anything for the exact eigendecomposition of rational matrices, so I doubt that anyone has written a program to do this before. However, you might have luck with the linear algebra features of sympy.


There is a "simple" method, simple in the sense that we need only few lines of procedure (here I use Maple, but Mathematica can do it also with PSLQ). Step 1: Maple calculates the roots of the characteristic polynomial with (for instance) 200 digits. Step 2: Maple uses the LLL method of Lenstra that gives the minimal polynomial of the roots ( If we know the roots with sufficient precision). Here the minimal polynomial has degree $1$ but, for sake of security, we ask for a polynomial of degree $2$. (if you obtain a polynomial of degree $2$, increase accuracy). The time of calculation, in the following random instance, is 0"12.

roll := rand(-10^10 .. 10^10): 
Digits:= 200:
f := 1:
 for i to 10 do f := f*(x-roll()*(1/roll())) end do;
f := numer(expand(f));
t := time():
 res := fsolve(f):

for i to 10 do 
print(MinimalPolynomial(res[i], 2)) end do;

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.