The shift-and-invert Arnoldi method for singular matrix pencils

Resumen Técnico: El Método de Arnoldi de Desplazamiento e Inversión para Paquetes de Matrices Singulares

Enunciado del Problema
El artículo aborda el desafío de resolver problemas de valores propios generalizados singulares, dispersos y de gran escala de la forma $Ax = \lambda Bx$, donde el paquete de matrices $A - \lambda B$ es singular (es decir, $\det(A - \lambda B) \equiv 0$). En tales casos, la definición clásica de valores propios falla porque cualquier escalar satisface la condición del determinante. Las soluciones significativas son los valores propios "regulares", definidos donde el rango de $A - \lambda_0 B$ desciende por debajo del rango normal del paquete.

Aunque existen métodos como el método escalera o la descomposición QZ para problemas pequeños y densos, son computacionalmente prohibitivos para sistemas dispersos de gran escala. Los enfoques recientes de Hochstenbach, Mehl y Plestenjak proponen regularizar paquetes singulares utilizando matrices densas elegidas aleatoriamente para separar los valores propios verdaderos de los espurios. Sin embargo, los autores argumentan que el uso de matrices de regularización densas destruye la estructura de dispersión del problema original, lo que conduce a sobrecargas significativas de rendimiento y memoria para aplicaciones de gran escala. Además, los solucionadores directos dispersos existentes no admiten nativamente la eliminación o adición de filas/columnas requerida para regularizar estos paquetes sin un desarrollo personalizado.

Metodología
Los autores proponen una estrategia determinista para regularizar paquetes singulares aprovechando la secuencia de pivoteo de una factorización LU dispersa. La metodología central implica los siguientes pasos:

1. Regularización mediante Pivoteo LU: En lugar de matrices densas aleatorias, el método utiliza una factorización LU dispersa de la matriz desplazada $A - \sigma B$ (donde $\sigma$ es un desplazamiento). Durante la factorización, si se encuentra un pivote cero (o cercano a cero), el algoritmo añade una fila y una columna a la matriz para mantener la factorización.

   • Paquete Aumentado (AGEP): El proceso construye matrices $F$ y $G$ tales que el paquete aumentado $\begin{bmatrix} A & G \\ F^* & 0 \end{bmatrix} - \lambda \begin{bmatrix} B & 0 \\ 0 & 0 \end{bmatrix}$ se vuelve regular. Las matrices $F$ y $G$ se derivan directamente de las matrices de permutación y del proceso de selección de pivotes de la factorización LU.
   • Paquete Proyectado (PGEP): Alternativamente, el método construye un paquete proyectado $G_\perp^* (A - \lambda B) F_\perp$ eliminando las filas y columnas correspondientes a los pivotes cero. $F_\perp$ y $G_\perp$ se derivan como la matriz identidad con filas específicas eliminadas.

2. Detección y Corrección de Rango: La factorización LU intenta naturalmente revelar el rango normal del paquete. Los autores señalan que, debido a errores de redondeo o configuraciones de tolerancia inadecuadas, el rango detectado ($\hat{k}$) podría diferir del rango normal verdadero ($k$).

   • Si el rango detectado es demasiado bajo (es decir, la tolerancia es demasiado agresiva), se propone un esquema de corrección de rango. Esto implica el uso de matrices ortogonales aleatorias para aumentar las $F$ y $G$ detectadas a la dimensión correcta, seguido de una actualización de bajo costo del complemento de Schur para resolver el sistema resultante.

3. Arnoldi de Desplazamiento e Inversión: El problema regularizado (ya sea AGEP o PGEP) se resuelve utilizando el método de Arnoldi de dos lados de desplazamiento e inversión. Los autores demuestran un teorema de equivalencia que muestra que la factorización de Arnoldi sobre el paquete aumentado (con un vector de inicio específico y un producto interno) produce los mismos valores de Ritz y la misma matriz de Hessenberg que el paquete proyectado. Esto permite el uso de la formulación proyectada para evitar la construcción explícita del sistema aumentado grande mientras se mantienen las mismas propiedades espectrales.

4. Problemas Rectangulares: El método se extiende a problemas de valores propios rectangulares ($A, B \in \mathbb{R}^{n \times m}$ con $n \neq m$). Para problemas rectangulares de rango completo, los solucionadores dispersos existentes con pivoteo de filas pueden utilizarse directamente para seleccionar un subproblema cuadrado y regular.

Contribuciones Clave

• Regularización Dispersa: El artículo introduce un método para generar matrices de regularización dispersas ($F, G, F_\perp, G_\perp$) directamente a partir de la secuencia de pivoteo de factorizaciones LU dispersas, preservando la dispersión del problema original.
• Detección de Rango: Demuestra que las factorizaciones LU dispersas estándar (utilizando pivoteo parcial, de alfil o completo) a menudo detectan con éxito el rango normal del paquete, facilitando la construcción automática de la regularización.
• Corrección de Rango: Se propone un algoritmo específico para corregir la regularización si la factorización LU subestima el rango normal, asegurando que el método permanezca robusto incluso con una detección de rango imperfecta.
• Equivalencia Teórica: El artículo establece un vínculo teórico riguroso entre los paquetes regularizados aumentados y proyectados en el contexto del método de Arnoldi, demostrando que producen aproximaciones espectrales idénticas bajo condiciones específicas.

Resultados
Se realizaron experimentos numéricos en tres tipos de problemas:

1. Actualización de Modelos de Elementos Finitos: Un modelo de estructura de celosía con 33 grados de libertad (resultando en un paquete de dimensión ~35,937). El método propuesto (utilizando pivoteo RLU) identificó correctamente el rango y calculó valores propios con residuos del orden de $10^{-14}$. En contraste, el enfoque aleatorio resultó en números de condición órdenes de magnitud más altos ($10^6$ vs $10^3$) y residuos más grandes.
2. Detección de Valores Propios Dobles: Un problema paramétrico que involucra la discretización de un Laplaciano 2D. El método basado en LU determinó correctamente el rango y logró alta precisión. El enfoque aleatorio requirió construir matrices densas, lo que llevó a costos de almacenamiento y computación significativamente más altos.
3. Problema de Valores Propios No Lineal Rectangular: Un problema rectangular de rango completo derivado de un problema de valores propios no lineal. El enfoque LU disperso fue significativamente más rápido (0.05s vs 2.9s para factorización y proyección) y requirió menos memoria (98MB vs 204MB) en comparación con el enfoque denso aleatorio.

En todos los casos, los métodos de regularización dispersa preservaron la dispersión de las matrices, lo que llevó a números de condición más bajos para los sistemas desplazados y una convergencia más rápida en comparación con la regularización densa aleatoria.

Significado y Afirmaciones
Los autores posicionan este trabajo como un concepto en lugar de una solución de software completa y optimizada, reconociendo que actualmente faltan solucionadores directos dispersos robustos para matrices singulares en las bibliotecas de software estándar. Sin embargo, el artículo afirma que:

• Rendimiento: Preservar la dispersión mediante pivoteo determinista es beneficioso tanto para el rendimiento del solucionador de valores propios como para la precisión de los resultados (mediante un mejor acondicionamiento).
• Viabilidad: Los solucionadores directos dispersos existentes pueden adaptarse para resolver estos problemas, particularmente para casos rectangulares de rango completo, sin necesidad de nueva infraestructura de software compleja.
• Robustez: El esquema de corrección de rango propuesto asegura que el método funcione incluso cuando la factorización LU no logra detectar perfectamente el rango normal.

El artículo concluye que, aunque el enfoque aleatorio es teóricamente sólido, el enfoque determinista de LU disperso ofrece una alternativa práctica, eficiente y precisa para problemas de valores propios singulares de gran escala, siempre que la brecha entre los valores singulares sea suficiente para la detección de rango. Se sugiere trabajo futuro para explorar solucionadores con pivoteo de volumen máximo global y para explotar estructuras de Kronecker en problemas multiparamétricos.