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. 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.