I am not sure if this is going to be faster than what you are doing now (which is already a clever method), but you can try the following.
- Compute $B = A^{-1}$. We know that $A$ has an eigenvalue in $(-1,1)$ iff $B$ has an eigenvalue outside $[-1,1].$
- You can run a few iterations of the power method to compute an approximate leading eigenvector of $B$: i.e., compute the first terms (let's say $m=10$ terms) of the sequence $v_0 = \text{random vector}$, $v_{k+1} = Bv_k$. The absolute value of the Rayleigh quotient $|v_m^* B v_m / v_m^* v_m|$ should converge to the spectral radius of $B$, for $m\to\infty$.
- If the resulting vector after $m$ iterations has a Rayleigh quotient $v_m^* B v_m / v_m^* v_m$ outside $[-1,1],$ then you can conclude for certainty that $B$ has an eigenvalue outside $[-1,1],$ since Rayleigh quotients lie in the convex hull of the eigenvalues of a matrix. This is easy to see by switching to an orthonormal basis in which $B$ is diagonal.
This gives you a quick test that may produce false negatives but not false positives, and thus it can be used to eliminate matrices quickly. Everything can be done in exact rational arithmetic.
Possible tricks for speedup (which might require testing):
- You can repeat the procedure for different random $v_0$ if you wish, but I imagine that it's more effective to increase the number of steps $m$ rather than restart the whole procedure.
- Instead of computing the inverse, you can compute a $LU$ or $LDL^*$ factorization of $A$, and use it to solve linear systems rather than computing inverses. This should be faster, in theory, but Sagemath is Python, so its performance is always tricky to predict.
- I am not expert in the area, but I think there are also LU-like factorization algorithms that can factor out denominators and work entirely in $\mathbb{Z}$. Using them might give further speedup.
- Alternatively, you can run all the steps in floating-point arithmetic, then replace $v_m$ with a nearby rational/integer vector $w_m$. If $w_m$, however computed, has $|w_m^* B w_m / w_m^*w_m| > 1$, then you can prove rigorously that $A$ has an eigenvector inside $(-1,1)$.
My guess is that this procedure is about as fast as computing a determinant (since they both require about $2/3 n^3$ operations asymptotically for an $n\times n$ matrix), but is more effective in weeding out matrices that have a small eigenvalue.
