Is there a fast way to check if a matrix has any small eigenvalues?

Question

I have hundreds of millions of symmetric 0/1-matrices of moderate size (say 20x20 to 30x30) which (obviously) have real eigenvalues.

I wish to extract from this list the tiny number of matrices that have no small eigenvalues, where I am calling an eigenvalue $\lambda$ small if it satisfies $|\lambda| < 1$.

Obviously I can do this by computing the characteristic polynomial of every matrix and proceeding from there. However I'd rather not calculate hundreds of millions of characteristic polynomials if I can avoid it.

So I am looking for techniques, heuristics etc that can be used to quickly test a matrix for the presence of a small eigenvalue so that I can immediately eliminate the matrix.

For example, if a matrix has determinant $0$ then it definitely has a small eigenvalue (namely $0$) and it can be eliminated immediately. Quick and dirty experiments with SageMath show that for a random sample of my matrices, the determinants can be calculated in about $6$% of the time required to compute the characteristic polynomials, and that about $15$% of the incoming matrices have determinant $0$. Therefore I gain a little bit by first computing determinants and throwing away those matrices with determinant $0$.

It seems though that this property is somewhat special because only a minuscule fraction of my hundreds of millions have no small eigenvalues. However I have searched fairly hard and have only found a few tangentially-related references to matrices with this property, though a couple more about matrices with no eigenvalues of modulus larger than one.

Question Are there are any other, perhaps more sophisticated, heuristics that can detect the presence of a small eigenvalue of a symmetric matrix?

Answer 1

Are you interested in a clever algorithm, or do you just want to get the answer fast?

If the latter, then I would suggest the following:

1. Use an established off-the-shelf eigenvalue solver.
2. Vectorize your code to work on many matrices at once.
3. Run it on a GPU.

There is an entire universe of experts (in both software and hardware) who dedicate their careers to creating libraries that do dense linear algebra operations extremely fast. Computing the eigenvalues for 100 million 30x30 matrices on the GPU should take less than one second. Computation time will be dominated by the cost of simply moving around the matrices in memory.

Here is some code using the python library jax which generates 1 million random symmetric 0-1 matrices in chunks of 100 thousand at a time, and checks which ones have small eigenvalues:

import numpy as np
import jax
import jax.numpy as jnp
from time import time
Generate 1 million random matrices in chunks of 100k
k = 30
chunk_size = int(1e5)
num_chunks = 10
t = time()
AA = []
for ii in range(num_chunks):
    A_nonsym = np.random.randint(0, 2, (chunk_size, k, k)) # random 0-1 matrices
    A = np.tril(A_nonsym) + np.tril(A_nonsym, -1).swapaxes(1,2) # symmetrize
    AA.append(A)
dt_generate_matrices = time() - t
print('dt_generate_matrices=', dt_generate_matrices)
Compute eigenvalues and check if small
bad_matrices_all_chunks = []
for ii, A in enumerate(AA):
    # move chunk of 100k matrices from memory to GPU
    t = time()
    A_gpu = jnp.array(A)
    dt_move = time() - t
# compute eigenvalues
t = time()
eigs = jnp.linalg.eigvalsh(A_gpu)
dt_eig = time() - t

# Check for small eigenvalues
t = time()
bad_matrices = jnp.any(np.abs(eigs) < 1.0, axis=1)
bad_matrices_all_chunks.append(bad_matrices)
dt_check = time() - t

print('chunk: ', ii, ', dt_move=', dt_move, ', dt_eig=', dt_eig)

When I run the above code, I get the following timing results:

dt_generate_matrices= 13.190782308578491
chunk:  0 , dt_move= 8.313929557800293 , dt_eig= 0.8072905540466309
chunk:  1 , dt_move= 0.2546844482421875 , dt_eig= 0.0008568763732910156
chunk:  2 , dt_move= 0.2372267246246338 , dt_eig= 0.0013964176177978516
chunk:  3 , dt_move= 0.23613929748535156 , dt_eig= 0.0001232624053955078
chunk:  4 , dt_move= 0.23468708992004395 , dt_eig= 0.00011324882507324219
chunk:  5 , dt_move= 0.23540091514587402 , dt_eig= 0.00011682510375976562
chunk:  6 , dt_move= 0.23517847061157227 , dt_eig= 0.00012135505676269531
chunk:  7 , dt_move= 0.23526215553283691 , dt_eig= 0.00012040138244628906
chunk:  8 , dt_move= 0.23635172843933105 , dt_eig= 0.00017189979553222656
chunk:  9 , dt_move= 0.23859858512878418 , dt_eig= 0.00011134147644042969

We see that generating 1 million random matrices took 13 seconds. Moving the first batch to the GPU and computing the eigenvalues took 9 seconds, which was almost all compiling time from jax's just-in-time compiler. After that, moving each chunk to the GPU took about 0.25 seconds. The actual eigenvalue computations for each chunk took a miniscule $10^{-4}$ seconds.

Multiplying by 100, the eigenvalue computation time for 100 million 30x30 matrices would be about 0.1 seconds (ignoring the time to generate the matrices, compile the code, and move data to/from the GPU).

Answer 2

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.

Answer 3

A symmetric matrix $A$ has an eigenvalue in $(-1,1)$ iff $A^2$ has an eigenvalue in $[0,1]$. Equivalently this means that $$ \min_{\Vert v\Vert=1} \frac{\Vert Av\Vert^2}{\Vert v\Vert^2} <1. $$ Hence you need to decide if the minimum of a nonnegative quadratic function on the unit sphere is $<1$.

I can think of two ways of doing this but I cannot comment on the computational complexity.

The first is a Monte Carlo type method. Sample $v$ uniformly on the sphere and compute $\Vert A v\Vert$. If the sample is large and there exists an eigenvalue of $A^2$ that is $<1$ you can pick detect this with high confidence. Technically, if the size of of the matrix is $n\times n$, then choose the coordinates $v_1,\dotsc, v_n$ of the random vector $v$ to be independent standard normal random variables and test if $$ \Vert A v\Vert^2 <\Vert v\Vert^2. $$

Here's a justifications. One can choose an orthonormal basis of your ambient space consisting of eigenvectors $e_1,\dotsc, e_n$ of $A^2$ corresponding to eigenvalues $\lambda_1\leq \cdots \leq \lambda_n$. If $\lambda_1<1$ then $\Vert Av\Vert^2<1$ in a neighborhood $U$ of $e_1$ on the unit sphere $\{\Vert v\Vert=1\}$. The probability that a random vector $v$ lands in $U$ is proportional to the size of $U$.
Here's a rough lower bound of the probability that $\Vert A\Vert<1$ given that $\Vert v\Vert=1$. Denote by $\sigma_k$ the area of the $k$-dimensional unit sphere.

Denote by $\lambda_\min$ respectively $\lambda_\max$ the smallest and respectively the largest of the eigenvalues of $A^2$. The nontrivial case is $\lambda_\min<1<\lambda_\max$, i.e., $\lambda_1=\lambda_\min$, $\lambda_n=\lambda_\max$.

Denote by $(v_1,\dotsc, v_n)$ the Euclidean coordinates determined by the basis $e_1,\dotsc, e_n$. Then $$ \Vert A v\Vert^2=\sum_{k=1}^n \lambda_kv_k^2. $$ if $\Vert v\Vert=1$, then $v_1^2=1-\sum_{k=2}^n v_k^2$ and $$ \Vert A v\Vert^2=\lambda_\min+\sum_{k=2}^n(\lambda_2-\lambda_\min)v_k^2 \leq\lambda_\min +(\lambda_\max-\lambda_\min)\sum_{k=2}v_k^2. $$ Thus, if $\Vert v\Vert =1$ and $$ \sum_{k=2}v_k^2<\frac{1-\lambda_\min}{\lambda_\max-\lambda_\min}=:\gamma(A), $$ then $\Vert Av\Vert<1$. In other words, if $\Vert v\Vert=1$ and $|v_1|>\sqrt{1-\gamma(A)}$, then $\Vert Av\Vert<1$. $\DeclareMathOperator{\Prob}{Prob}$. Hence $$ \Prob\big[ \Vert Av\Vert<1\big|\;\Vert v\Vert=1\big]\geq \frac{1}{\sigma_{n-1}}{\rm area}\;\big(\{v\in S^{n-1};\;v_1^2>1-\gamma(A)\}\big) $$ $$ =\frac{\sigma_{n-2}}{\sigma_{n-1}}\int_{1-\gamma(A)\leq v_1^2\leq 1}^1 (1-v_1^2)^{(n-3)/2} dv_1=\frac{\sigma_{n-2}}{\sigma_{n-1}}\int_{1-\gamma(A)}^1 (1-s)^{(n-3)/2}s^{-1/2} ds. $$ (The component $v_1$, as a random variable defined on the unit sphere has a Beta distribution) $$ \frac{\sigma_{n-2}}{\sigma_{n-1}}\int_{1-\gamma(A)}^1 (1-s)^{(n-3)/2}s^{-1/2} ds=\frac{1}{B((n-1)/2,1/2)}\int_0^{\gamma(A)} s^{(n-3)/2}(1-s)^{-1/2} ds. $$

We have $$ B((n-1)/2,1/2)\sim \Gamma(1/2)\sqrt{\frac{2}{n-1}}, $$ and $$ \int_0^{\gamma(A)} s^{(n-3)/2}(1-s)^{-1/2} ds \approx \frac{2}{n-1}\gamma(A)^{\frac{n-1}{2}}.$$

The expected time to detect $\Vert Av\Vert<1$ is about $\frac{1}{\sqrt{n}}\gamma(A)^{-n/2}$

Hence for $n<30$ the sample size ought to be in the millions. I cannot estimate how long it would take to run such samples. (I've run samples of this size on my laptop.) The computation of $\Vert Av\Vert$ is not very time consuming since $A$ is a $0/1$ matrix.

If one your samples yield a vector such that $\Vert Av\Vert^2<1$ then game over.

If your samples yield only numbers $\Vert Av\Vert^2$ substantially bigger than $1$ then you can say with high confidence that $A$ has no small eigenvalues. If the numbers $\Vert A v\Vert$ are bigger than $1$ "many" are close to $1$ there is a bit of ambiguity.

Question. A natural question comes to mind. Suppose that $$ A=(A_{ij})_{1\leq i,j\leq n} $$ is a random symmetric $0/1$ matrix: the coefficients $A_{ij}$, $1\leq i\leq j\leq n$ are independent Bernoulli random variables with probability of success $1/2$. It would be interesting to know the probaility $p(n)$ that $A$ has an eigenvalue in $[-1,1]$

The second approach is by gradient descent. Consider the function $$ f:\{\Vert v\Vert=1\}\to [0,\infty),\;\;f(v)=\Vert Av\Vert^2. $$ A flow line $v(t)$ of the negative gradient flow $$ \frac{dv}{dt}=-\nabla f(v) $$ will converge exponentially to an eigenvector of $A^2$. If the initial condition $v(0)$ is uniformly random on the unit sphere, then, with probability $1$, this flow line will converge exponentially to an eigenvector corresponding to a minimal eigenvalue.

Solve numerically this equation with random initial condition $v(0)$. I speculate that this discretization will lead you very fast to a decision concerning small e-values. The speed of convergence depends on how "packed" are the eigenvalues of $A^2$: if they are all packed in a small interval, the convergence will be slower.

Again, I cannot comment on the computational complexity of this approach.

Answer 4

This subtle question requires subtle answers. To begin with, let me discard a few strategies that seem erroneous when the matrices $A\in{\bf Sym}_n$ have a very large size $n$.

• Random sampling: compute $x^TAx$ with $x$ randomly chosen in the the unit sphere. It was shown (in T. Gallay & D. S., Comm. Pure Appl. Math. 2012) that this number is, with very high probability, close to $\frac1n{\rm Tr}A$, and thus does not give information about the small eigenvalues.
• Solve the differential equation $\dot x=-Ax$ (assuming that $A$ is positive). Doing so with precision requires the inversion of some matrix $hI_n+A$, which is way too costly.

Is suggest instead the following, based upon the QR method (see my book Matrices, 2nd edition, chapter 13.2). In principle, QR is meant to compute the whole spectrum of $A$, which is not needed here, but we only have to implement a small number of iterations (see below). One of the advantages of QR is that the matrix $A$ is transformed by orthogonal conjugations, and thus the numerical errors are not amplified: QR is a very stable method.

First of all, you may assume that $A$ is positive (otherwise replace $A$ by $A^2$). Then put $A$ in tridiagonal form, in about $\frac23n^3$ operations. Every QR iteration will cost only $O(n)$ iterations and the iterates $A^{(k)}$ remain tridiagonal ; this is incredibly cheap ! Then proceed: it is known that the QR sequence converges to ${\rm diag}(\mu_1,\ldots,\mu_n)$ where $\mu_1\ge\cdots\ge\mu_n(\ge0)$ are the eigenvalues. You are interested in $\mu_n$, aren't you ? Good news: the convergence is the fastest (by far) at the south-east corner. In other words you need very few iterations (say a dozen) to have a very good approximation of $\mu_n$.

Answer 5

You can check if the matrix $$A^T A - I$$ is positive definite. This can be done with Cholesky Decomposition, which runs in time $O(n^3)$ and fails for non-positive-definite matrices.

If $A^TA -I$ is not positive definite, then there is some $x\in\mathbb{R}^n$ such that $x(A^T A -I )x < 0$, therefore $$x^T(A^T A)x< \Vert x\Vert^2 \implies \lambda_\min(A^TA)<1\implies \lambda_\min(A)\in(-1,1)$$

---

I came across something surprising: by generating random 0/1 matrices, it seems that all of them contain an eigenvalue $\lambda_\min(A) \le [-1,1]$. Is it just that counter-examples are rare?

Answer 6

YES, the pari/gp function called qfsign() does use rational numbers only!!!!! You may use it and not fear!

Got a reply from Bill Allombert himself, on a mailing list

Pari-gp has a version of what I did, naming it after Gauss rather than Lagrange, in the French manner. It is quick, just field operations, so all are rational numbers if you give it integers.

https://fr.wikipedia.org/wiki/R%C3%A9duction_de_Gauss

because I can then copy and paste each line and not introduce errors.

Something I have programmed: see comment of Brendan McKay, but as an if and only if... Sylvester's Law of Inertia says that we may decide this if the counts of negative, positive, zero eigenvalues change a little among three symmetric matrices, call them $H, H+I, H-I.$

The hard part was making it give Latex code as output. This is all C++ with GMP; note I wrote my own fractions ... Altogether it is three files, about 200 lines each: myfrac.h /// matrixrational.h /// matrix_congruence.cc

Next given an integer symmetric matrix $H,$ we may solve $P^T H P = D$ diagonal for $P$ nonsingular with rational entries. Once we have $D$ we just count the diagonal elements, how many positive, how many negative, how many zero. The same counts doing the same thing for $H-I$ and $H+I$

$$ H = \left( \begin{array}{rrrrr} 2 & - 1 & - 1 & - 1 & - 1 \\ - 1 & 2 & 0 & 0 & 0 \\ - 1 & 0 & 2 & 0 & 0 \\ - 1 & 0 & 0 & 2 & 0 \\ - 1 & 0 & 0 & 0 & 2 \\ \end{array} \right) $$ $$ D_0 = H $$ $$ E_j^T D_{j-1} E_j = D_j $$ $$ P_{j-1} E_j = P_j $$ $$ E_j^{-1} Q_{j-1} = Q_j $$ $$ P_j Q_j = Q_j P_j = I $$ $$ P_j^T H P_j = D_j $$ $$ Q_j^T D_j Q_j = H $$

$$ H = \left( \begin{array}{rrrrr} 2 & - 1 & - 1 & - 1 & - 1 \\ - 1 & 2 & 0 & 0 & 0 \\ - 1 & 0 & 2 & 0 & 0 \\ - 1 & 0 & 0 & 2 & 0 \\ - 1 & 0 & 0 & 0 & 2 \\ \end{array} \right) $$

$$ \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc $$

$$ P^T H P = D $$

$$\left( \begin{array}{rrrrr} 1 & 0 & 0 & 0 & 0 \\ \frac{ 1 }{ 2 } & 1 & 0 & 0 & 0 \\ \frac{ 2 }{ 3 } & \frac{ 1 }{ 3 } & 1 & 0 & 0 \\ 1 & \frac{ 1 }{ 2 } & \frac{ 1 }{ 2 } & 1 & 0 \\ 2 & 1 & 1 & 1 & 1 \\ \end{array} \right) \left( \begin{array}{rrrrr} 2 & - 1 & - 1 & - 1 & - 1 \\ - 1 & 2 & 0 & 0 & 0 \\ - 1 & 0 & 2 & 0 & 0 \\ - 1 & 0 & 0 & 2 & 0 \\ - 1 & 0 & 0 & 0 & 2 \\ \end{array} \right) \left( \begin{array}{rrrrr} 1 & \frac{ 1 }{ 2 } & \frac{ 2 }{ 3 } & 1 & 2 \\ 0 & 1 & \frac{ 1 }{ 3 } & \frac{ 1 }{ 2 } & 1 \\ 0 & 0 & 1 & \frac{ 1 }{ 2 } & 1 \\ 0 & 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 0 & 1 \\ \end{array} \right) = \left( \begin{array}{rrrrr} 2 & 0 & 0 & 0 & 0 \\ 0 & \frac{ 3 }{ 2 } & 0 & 0 & 0 \\ 0 & 0 & \frac{ 4 }{ 3 } & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ \end{array} \right) $$ $$ $$

$$ Q^T D Q = H $$

$$\left( \begin{array}{rrrrr} 1 & 0 & 0 & 0 & 0 \\ - \frac{ 1 }{ 2 } & 1 & 0 & 0 & 0 \\ - \frac{ 1 }{ 2 } & - \frac{ 1 }{ 3 } & 1 & 0 & 0 \\ - \frac{ 1 }{ 2 } & - \frac{ 1 }{ 3 } & - \frac{ 1 }{ 2 } & 1 & 0 \\ - \frac{ 1 }{ 2 } & - \frac{ 1 }{ 3 } & - \frac{ 1 }{ 2 } & - 1 & 1 \\ \end{array} \right) \left( \begin{array}{rrrrr} 2 & 0 & 0 & 0 & 0 \\ 0 & \frac{ 3 }{ 2 } & 0 & 0 & 0 \\ 0 & 0 & \frac{ 4 }{ 3 } & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ \end{array} \right) \left( \begin{array}{rrrrr} 1 & - \frac{ 1 }{ 2 } & - \frac{ 1 }{ 2 } & - \frac{ 1 }{ 2 } & - \frac{ 1 }{ 2 } \\ 0 & 1 & - \frac{ 1 }{ 3 } & - \frac{ 1 }{ 3 } & - \frac{ 1 }{ 3 } \\ 0 & 0 & 1 & - \frac{ 1 }{ 2 } & - \frac{ 1 }{ 2 } \\ 0 & 0 & 0 & 1 & - 1 \\ 0 & 0 & 0 & 0 & 1 \\ \end{array} \right) = \left( \begin{array}{rrrrr} 2 & - 1 & - 1 & - 1 & - 1 \\ - 1 & 2 & 0 & 0 & 0 \\ - 1 & 0 & 2 & 0 & 0 \\ - 1 & 0 & 0 & 2 & 0 \\ - 1 & 0 & 0 & 0 & 2 \\ \end{array} \right) $$

$$ \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc $$

next $H-I$ $$ P^T H P = D $$

$$\left( \begin{array}{rrrrr} 1 & 0 & 0 & 0 & 0 \\ 2 & 1 & 1 & 0 & 0 \\ 0 & - \frac{ 1 }{ 2 } & \frac{ 1 }{ 2 } & 0 & 0 \\ - 1 & - 1 & - 1 & 1 & 0 \\ - \frac{ 1 }{ 2 } & - \frac{ 1 }{ 2 } & - \frac{ 1 }{ 2 } & - \frac{ 1 }{ 2 } & 1 \\ \end{array} \right) \left( \begin{array}{rrrrr} 1 & - 1 & - 1 & - 1 & - 1 \\ - 1 & 1 & 0 & 0 & 0 \\ - 1 & 0 & 1 & 0 & 0 \\ - 1 & 0 & 0 & 1 & 0 \\ - 1 & 0 & 0 & 0 & 1 \\ \end{array} \right) \left( \begin{array}{rrrrr} 1 & 2 & 0 & - 1 & - \frac{ 1 }{ 2 } \\ 0 & 1 & - \frac{ 1 }{ 2 } & - 1 & - \frac{ 1 }{ 2 } \\ 0 & 1 & \frac{ 1 }{ 2 } & - 1 & - \frac{ 1 }{ 2 } \\ 0 & 0 & 0 & 1 & - \frac{ 1 }{ 2 } \\ 0 & 0 & 0 & 0 & 1 \\ \end{array} \right) = \left( \begin{array}{rrrrr} 1 & 0 & 0 & 0 & 0 \\ 0 & - 2 & 0 & 0 & 0 \\ 0 & 0 & \frac{ 1 }{ 2 } & 0 & 0 \\ 0 & 0 & 0 & 2 & 0 \\ 0 & 0 & 0 & 0 & \frac{ 3 }{ 2 } \\ \end{array} \right) $$ $$ $$

$$ Q^T D Q = H $$

$$\left( \begin{array}{rrrrr} 1 & 0 & 0 & 0 & 0 \\ - 1 & \frac{ 1 }{ 2 } & - 1 & 0 & 0 \\ - 1 & \frac{ 1 }{ 2 } & 1 & 0 & 0 \\ - 1 & 1 & 0 & 1 & 0 \\ - 1 & 1 & 0 & \frac{ 1 }{ 2 } & 1 \\ \end{array} \right) \left( \begin{array}{rrrrr} 1 & 0 & 0 & 0 & 0 \\ 0 & - 2 & 0 & 0 & 0 \\ 0 & 0 & \frac{ 1 }{ 2 } & 0 & 0 \\ 0 & 0 & 0 & 2 & 0 \\ 0 & 0 & 0 & 0 & \frac{ 3 }{ 2 } \\ \end{array} \right) \left( \begin{array}{rrrrr} 1 & - 1 & - 1 & - 1 & - 1 \\ 0 & \frac{ 1 }{ 2 } & \frac{ 1 }{ 2 } & 1 & 1 \\ 0 & - 1 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & \frac{ 1 }{ 2 } \\ 0 & 0 & 0 & 0 & 1 \\ \end{array} \right) = \left( \begin{array}{rrrrr} 1 & - 1 & - 1 & - 1 & - 1 \\ - 1 & 1 & 0 & 0 & 0 \\ - 1 & 0 & 1 & 0 & 0 \\ - 1 & 0 & 0 & 1 & 0 \\ - 1 & 0 & 0 & 0 & 1 \\ \end{array} \right) $$

$$ \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc $$

next $H+I$

$$ P^T H P = D $$

$$\left( \begin{array}{rrrrr} 1 & 0 & 0 & 0 & 0 \\ \frac{ 1 }{ 3 } & 1 & 0 & 0 & 0 \\ \frac{ 3 }{ 8 } & \frac{ 1 }{ 8 } & 1 & 0 & 0 \\ \frac{ 3 }{ 7 } & \frac{ 1 }{ 7 } & \frac{ 1 }{ 7 } & 1 & 0 \\ \frac{ 1 }{ 2 } & \frac{ 1 }{ 6 } & \frac{ 1 }{ 6 } & \frac{ 1 }{ 6 } & 1 \\ \end{array} \right) \left( \begin{array}{rrrrr} 3 & - 1 & - 1 & - 1 & - 1 \\ - 1 & 3 & 0 & 0 & 0 \\ - 1 & 0 & 3 & 0 & 0 \\ - 1 & 0 & 0 & 3 & 0 \\ - 1 & 0 & 0 & 0 & 3 \\ \end{array} \right) \left( \begin{array}{rrrrr} 1 & \frac{ 1 }{ 3 } & \frac{ 3 }{ 8 } & \frac{ 3 }{ 7 } & \frac{ 1 }{ 2 } \\ 0 & 1 & \frac{ 1 }{ 8 } & \frac{ 1 }{ 7 } & \frac{ 1 }{ 6 } \\ 0 & 0 & 1 & \frac{ 1 }{ 7 } & \frac{ 1 }{ 6 } \\ 0 & 0 & 0 & 1 & \frac{ 1 }{ 6 } \\ 0 & 0 & 0 & 0 & 1 \\ \end{array} \right) = \left( \begin{array}{rrrrr} 3 & 0 & 0 & 0 & 0 \\ 0 & \frac{ 8 }{ 3 } & 0 & 0 & 0 \\ 0 & 0 & \frac{ 21 }{ 8 } & 0 & 0 \\ 0 & 0 & 0 & \frac{ 18 }{ 7 } & 0 \\ 0 & 0 & 0 & 0 & \frac{ 5 }{ 2 } \\ \end{array} \right) $$ $$ $$

$$ Q^T D Q = H $$

$$\left( \begin{array}{rrrrr} 1 & 0 & 0 & 0 & 0 \\ - \frac{ 1 }{ 3 } & 1 & 0 & 0 & 0 \\ - \frac{ 1 }{ 3 } & - \frac{ 1 }{ 8 } & 1 & 0 & 0 \\ - \frac{ 1 }{ 3 } & - \frac{ 1 }{ 8 } & - \frac{ 1 }{ 7 } & 1 & 0 \\ - \frac{ 1 }{ 3 } & - \frac{ 1 }{ 8 } & - \frac{ 1 }{ 7 } & - \frac{ 1 }{ 6 } & 1 \\ \end{array} \right) \left( \begin{array}{rrrrr} 3 & 0 & 0 & 0 & 0 \\ 0 & \frac{ 8 }{ 3 } & 0 & 0 & 0 \\ 0 & 0 & \frac{ 21 }{ 8 } & 0 & 0 \\ 0 & 0 & 0 & \frac{ 18 }{ 7 } & 0 \\ 0 & 0 & 0 & 0 & \frac{ 5 }{ 2 } \\ \end{array} \right) \left( \begin{array}{rrrrr} 1 & - \frac{ 1 }{ 3 } & - \frac{ 1 }{ 3 } & - \frac{ 1 }{ 3 } & - \frac{ 1 }{ 3 } \\ 0 & 1 & - \frac{ 1 }{ 8 } & - \frac{ 1 }{ 8 } & - \frac{ 1 }{ 8 } \\ 0 & 0 & 1 & - \frac{ 1 }{ 7 } & - \frac{ 1 }{ 7 } \\ 0 & 0 & 0 & 1 & - \frac{ 1 }{ 6 } \\ 0 & 0 & 0 & 0 & 1 \\ \end{array} \right) = \left( \begin{array}{rrrrr} 3 & - 1 & - 1 & - 1 & - 1 \\ - 1 & 3 & 0 & 0 & 0 \\ - 1 & 0 & 3 & 0 & 0 \\ - 1 & 0 & 0 & 3 & 0 \\ - 1 & 0 & 0 & 0 & 3 \\ \end{array} \right) $$

$$ \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc $$

Answer 7

I am not currently allowed to comment, so I am writing this trivial observation as an answer instead.

A very slight variation of the method you are currently using is possible. The determinant is the product of the eigenvalues and it is easy to see that a product of real numbers can only be in the interval $(-1,1)$ if at least one of the numbers is also in $(-1,1)$. Therefore rather than rejecting matrices with determinant equal to zero you can reject all those with determinant in $(-1,1)$.

Answer 8

The power iteration technique is often used to calculate the dominant eigenvalue/eigenvector of an arbitrary matrix. The power iteration technique is simple and efficient, but it does not always work and it only gives the dominant eigenvalue/eigenvector, but power iteration works better for symmetric matrices. We can use the power iteration (or some other iterative algorithm) to calculate the minimum eigenvalue of a symmetric matrix as well.

Suppose that $A$ is symmetric (the eigenvalues need to be real for this to work). Let $B=A^2$. Suppose that $B$ has eigenvalues $\lambda_1\geq\dots\geq\lambda_n$. Let $\mu>0$. Let $C=\mu\cdot I-B$. Then $C$ has eigenvalues $\mu-\lambda_1,\dots,\mu-\lambda_n$. If $|\mu-\lambda_n|\geq|\mu-\lambda_1|$, then the spectral radius of $C$ will be $\mu-\lambda_n$ and this spectral radius may be computed using power iteration. While we want $|\mu-\lambda_n|\geq|\mu-\lambda_1|$, it is better if $\mu$ is not too large since we want the power iteration to converge quickly.

To use the power iteration, we begin with a random vector $v_0$, and then we set $v_{r+1}=Cv_r$ for all $r$. Each step in the power iteration can be computed efficiently since $Cv_r=\mu\cdot v_r-A(Av_r)$ where the parentheses indicate the fastest order of operations (we do not need to compute $B$ to get $Cv_r$). In this case, $\|v_{r+1}\|/\|v_r\|$ will converge to the spectral radius of $C$.