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 '13 at 18:20

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

