¿Existe una forma rápida de verificar si una matriz tiene valores propios pequeños?

Question

Tengo cientos de millones de matrices simétricas de 0/1 de tamaño moderado (digamos 20x20 a 30x30) que (obviamente) tienen eigenvalores reales.

Deseo extraer de esta lista el diminuto número de matrices que tienen ningún eigenvalor pequeño, donde estoy llamando eigenvalor $\lambda$ pequeño si satisface $|\lambda| < 1$.

Obviamente puedo hacer esto calculando el polinomio característico de cada matriz y procediendo a partir de ahí. Sin embargo, preferiría no calcular cientos de millones de polinomios característicos si puedo evitarlo.

Por lo tanto, estoy buscando técnicas, heurísticas, etc. que se pueden utilizar para probar rápidamente una matriz en busca de la presencia de un eigenvalor pequeño para poder eliminar inmediatamente la matriz.

Por ejemplo, si una matriz tiene determinante $0$ entonces definitivamente tiene un eigenvalor pequeño (o sea, $0$) y se puede eliminar inmediatamente. Experimentos rápidos y sucios con SageMath muestran que para una muestra aleatoria de mis matrices, los determinantes se pueden calcular en aproximadamente $6$% del tiempo requerido para calcular los polinomios característicos, y que alrededor del $15$% de las matrices entrantes tienen determinante $0$. Por lo tanto, gano un poco al calcular primero los determinantes y desechar esas matrices con determinante $0$.

Parece que esta propiedad es algo especial porque solo una fracción minúscula de mis cientos de millones no tienen eigenvalores pequeños. Sin embargo, he buscado bastante y solo he encontrado algunas referencias tangencialmente relacionadas con matrices con esta propiedad, aunque un par más sobre matrices sin eigenvalores de módulo mayor que uno.

Pregunta ¿Hay alguna otra heurística, tal vez más sofisticada, que pueda detectar la presencia de un eigenvalor pequeño de una matriz simétrica?

Answer 1

¿Estás interesado en un algoritmo ingenioso, o simplemente quieres obtener la respuesta rápido?

Si es lo último, entonces te sugeriría lo siguiente:

1. Usa un solucionador establecido de valores propios.
2. Vectoriza tu código para trabajar en muchas matrices a la vez.
3. Ejecútalo en una GPU.

Existe un universo entero de expertos (tanto en software como en hardware) que dedican sus carreras a crear bibliotecas que realizan operaciones de álgebra lineal densa extremadamente rápido. Calcular los valores propios de 100 millones de matrices de 30x30 en la GPU debería tomar menos de un segundo. El tiempo de computación estará dominado por el costo de simplemente mover las matrices en la memoria.

Aquí tienes un código que utiliza la librería de Python jax que genera 1 millón de matrices aleatorias simétricas de 0 a 1 en grupos de 100 mil a la vez, y verifica cuáles tienen valores propios pequeños:

import numpy as np
import jax
import jax.numpy as jnp
from time import time

#### Generar 1 millón de matrices aleatorias en grupos de 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)) # matrices aleatorias de 0 a 1
    A = np.tril(A_nonsym) + np.tril(A_nonsym, -1).swapaxes(1,2) # simetrizar
    AA.append(A)
dt_generate_matrices = time() - t
print('dt_generate_matrices=', dt_generate_matrices)

#### Calcular valores propios y verificar si son pequeños ####
bad_matrices_all_chunks = []
for ii, A in enumerate(AA):
    # mover grupo de 100k matrices de la memoria a la GPU
    t = time()
    A_gpu = jnp.array(A)
    dt_move = time() - t

    # calcular valores propios
    t = time()
    eigs = jnp.linalg.eigvalsh(A_gpu)
    dt_eig = time() - t

    # Verificar valores propios pequeños
    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('grupo: ', ii, ', dt_move=', dt_move, ', dt_eig=', dt_eig)

Cuando ejecuto el código anterior, obtengo los siguientes resultados de tiempo:

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

Vemos que generar 1 millón de matrices aleatorias tomó 13 segundos. Mover el primer grupo a la GPU y calcular los valores propios tomó 9 segundos, que fue casi todo tiempo de compilación del compilador en tiempo real de jax. Después de eso, mover cada grupo a la GPU tomó alrededor de 0,25 segundos. Los cálculos reales de los valores propios para cada grupo tomaron un tiempo ínfimo de $10^{-4}$ segundos.

Multiplicando por 100, el tiempo de cálculo de valores propios para 100 millones de matrices de 30x30 sería de alrededor de 0,1 segundos (ignorando el tiempo para generar las matrices, compilar el código y mover datos hacia/desde la GPU).

Answer 2

No estoy seguro de si esto será más rápido que lo que estás haciendo ahora (que ya es un método inteligente), pero puedes intentarlo.

• Calcular $B = A^{-1}$. Sabemos que $A$ tiene un valor propio en $(-1,1)$ si y solo si $B$ tiene un valor propio fuera de $[-1,1].$
• Puedes ejecutar algunas iteraciones del método de la potencia para calcular un vector propio principal aproximado de $B$: es decir, calcular los primeros términos (digamos $m=10$ términos) de la secuencia $v_0 = \text{vector aleatorio}$, $v_{k+1} = Bv_k$. El valor absoluto del cociente de Rayleigh $|v_m^* B v_m / v_m^* v_m|$ debería converger al radio espectral de $B$, para $m\to\infty$.
• Si el vector resultante después de $m$ iteraciones tiene un cociente de Rayleigh $v_m^* B v_m / v_m^* v_m$ fuera de $[-1,1],$ entonces puedes concluir con certeza que $B$ tiene un valor propio fuera de $[-1,1],$ ya que los cocientes de Rayleigh están en el envoltorio convexo de los valores propios de una matriz. Esto es fácil de ver al cambiar a una base ortonormal en la que $B$ es diagonal.

Esto te da una prueba rápida que puede producir falsos negativos pero no falsos positivos, y por lo tanto puede usarse para eliminar matrices rápidamente. Todo se puede hacer en aritmética racional exacta.

Posibles trucos para acelerar (que pueden requerir pruebas):

• Puedes repetir el procedimiento para diferentes $v_0$ aleatorios si quieres, pero imagino que es más efectivo aumentar el número de pasos $m$ en lugar de reiniciar todo el procedimiento.
• En lugar de calcular la inversa, puedes calcular una factorización $LU$ o $LDL^*$ de $A$, y usarla para resolver sistemas lineales en lugar de calcular inversas. Esto debería ser más rápido, en teoría, pero Sagemath es Python, así que su rendimiento siempre es complicado de predecir.
• No soy experto en el área, pero creo que también hay algoritmos de factorización tipo LU que pueden eliminar denominadores y trabajar totalmente en $\mathbb{Z}$. Usarlos podría dar un mayor aumento de velocidad.
• Alternativamente, puedes ejecutar todos los pasos en aritmética de punto flotante, luego reemplazar $v_m$ con un vector racional/enterto cercano $w_m$. Si $w_m$, calculado de cualquier manera, tiene $|w_m^* B w_m / w_m^*w_m| > 1$, entonces puedes demostrar rigurosamente que $A$ tiene un vector propio dentro de $(-1,1)$.

Supongo que este procedimiento es aproximadamente tan rápido como calcular un determinante (ya que ambos requieren aproximadamente $2/3 n^3$ operaciones asintóticamente para una matriz $n\times n$), pero es más efectivo para filtrar matrices que tienen un valor propio pequeño.

Answer 3

Esta sutil pregunta requiere respuestas sutiles. Para empezar, descartemos algunas estrategias que parecen erróneas cuando las matrices $A\in{\bf Sym}_n$ tienen un tamaño muy grande $n$.

• Muestreo aleatorio: calcular $x^TAx$ con $x$ elegido al azar en la esfera unitaria. Se demostró (en T. Gallay & D. S., Comm. Pure Appl. Math. 2012) que este número es, con una probabilidad muy alta, cercano a $\frac1n{\rm Tr}A$, y por lo tanto no proporciona información sobre los pequeños autovalores.
• Resolver la ecuación diferencial $\dot x=-Ax$ (asumiendo que $A$ es positiva). Hacerlo con precisión requiere la inversión de alguna matriz $hI_n+A$, lo cual es demasiado costoso.

Sugiero en su lugar lo siguiente, basado en el método QR (ver mi libro Matrices, 2da edición, capítulo 13.2). En principio, QR está destinado a calcular todo el espectro de $A$, lo cual no es necesario aquí, pero solo tenemos que implementar un pequeño número de iteraciones (ver abajo). Una de las ventajas de QR es que la matriz $A$ es transformada por conjugaciones ortogonales, y así los errores numéricos no se amplifican: QR es un método muy estable.

En primer lugar, puedes suponer que $A$ es positiva (de lo contrario reemplaza $A$ por $A^2$). Luego pon $A$ en forma tridiagonal, en aproximadamente $\frac23n^3$ operaciones. Cada iteración de QR solo costará $O(n)$ iteraciones y los iterados $A^{(k)}$ seguirán siendo tridiagonales; ¡esto es increíblemente barato! Luego continúa: se sabe que la secuencia QR converge a ${\rm diag}(\mu_1,\ldots,\mu_n)$ donde $\mu_1\ge\cdots\ge\mu_n(\ge0)$ son los autovalores. ¿Estás interesado en $\mu_n$? Buenas noticias: la convergencia es la más rápida (por mucho) en la esquina sureste. En otras palabras, necesitas muy pocas iteraciones (digamos una docena) para tener una muy buena aproximación de $\mu_n$.

Answer 4

Puedes verificar si la matriz $$A^T A - I$$ es definitiva positiva. Esto se puede hacer con la Descomposición de Cholesky, que se ejecuta en tiempo $O(n^3)$ y falla para matrices no definitivas positivas.

Si $A^TA -I$ no es definitiva positiva, entonces existe algún $x\in\mathbb{R}^n$ tal que $x(A^T A -I )x < 0$, por lo tanto $$x^T(A^T A)x< \Vert x\Vert^2 \implies \lambda_\min(A^TA)<1\implies \lambda_\min(A)\in(-1,1)$$

---

Me encontré con algo sorprendente: al generar matrices aleatorias de 0/1, parece que todas contienen un valor propio $\lambda_\min(A) \le [-1,1]$. ¿Es solo que los contraejemplos son raros?

Answer 5

Actualmente no se me permite comentar, así que escribo esta observación trivial como respuesta en su lugar.

Es posible una variación muy ligera del método que estás usando actualmente. El determinante es el producto de los autovalores y es fácil ver que un producto de números reales solo puede estar en el intervalo $(-1,1)$ si al menos uno de los números también está en $(-1,1)$. Por lo tanto, en lugar de rechazar las matrices con determinante igual a cero, puedes rechazar todas aquellas con un determinante en $(-1,1)$.