The eigenvalue decomposition of normal matrices by the skew-symmetric part

1. Planteamiento del Problema

El cálculo de la descomposición en valores propios (EVD) y la descomposición de Schur real (RSD) es un problema fundamental en álgebra lineal. Para matrices normales reales ($AA^\top = A^\top A$), se sabe que son diagonalizables mediante transformaciones de similitud unitarias (ortogonales en el caso real), lo que implica que su forma de Schur es diagonal (o bloques diagonales de $2\times2$ para valores propios complejos conjugados).

Sin embargo, los algoritmos estándar para matrices generales (como el algoritmo QR con reducción a forma de Hessenberg) no aprovechan la estructura especial de las matrices normales, operando sobre matrices densas y con un costo computacional elevado ($O(n^3)$) sin garantizar la diagonalización final de manera eficiente.

El problema específico abordado es: ¿Cómo calcular la EVD/RSD de una matriz normal real densa de manera más eficiente, especialmente cuando la matriz tiene pocos valores propios reales (como es común en matrices ortogonales aleatorias)?

2. Metodología Propuesta

Los autores proponen un nuevo método que explota la estructura de la parte antisimétrica de la matriz normal, denotada como $\Omega = \text{skew}(A) = \frac{1}{2}(A - A^\top)$.

Fundamentos Teóricos Clave

1. Conmutatividad: Para una matriz normal $A$, su parte simétrica $H$ y su parte antisimétrica $\Omega$ conmutan ($H\Omega = \Omega H$). Esto implica que son diagonalizables simultáneamente.
2. Relación de Subespacios Invariantes:
   • Los subespacios invariantes de $\Omega$ coinciden con los subespacios invariantes de $A$ para los valores propios con partes imaginarias no nulas y distintas.
   • Si $A$ tiene $r$ valores propios reales, el núcleo de $\Omega$ corresponde exactamente al subespacio de los valores propios reales de $A$.
   • Si hay valores propios complejos con partes imaginarias repetidas (agrupadas), se requiere un paso adicional de refinamiento.

El Algoritmo (Aritmética Exacta y de Punto Flotante)

El método se estructura en los siguientes pasos principales (Algoritmo 4.1 en el papel):

1. Cálculo de la Parte Antisimétrica: Se calcula $\Omega = \text{skew}(A)$.
2. Reducción Tridiagonal: Se reduce $\Omega$ a forma tridiagonal antisimétrica utilizando reflectores de Householder.
3. Permutación y SVD: Se aplica una permutación "par-impar" (even-odd) para transformar la matriz tridiagonal en una forma que revela un problema de valores singulares bidiagonal de tamaño reducido ($\approx n/2$). Se calcula la SVD de esta matriz bidiagonal.
4. Identificación de Clústeres: Se analizan los valores singulares de $\Omega$ (que corresponden a las partes imaginarias de los valores propios de $A$):
   • Caso 1 (Partes imaginarias distintas): Los vectores de Schur de $\Omega$ son directamente los vectores de Schur de $A$. No se requiere más trabajo.
   • Caso 2 (Partes imaginarias repetidas/clústeres): Si hay multiplicidad en las partes imaginarias, se extrae el subespacio correspondiente, se proyecta la matriz original $A$ sobre este subespacio y se resuelve un problema de valores propios simétrico-skew-Hamiltoniano de menor dimensión.
   • Caso 3 (Valores propios reales): Si $\Omega$ tiene valores singulares cero (núcleo no trivial), esto corresponde a los valores propios reales de $A$. Se proyecta $A$ sobre este núcleo y se resuelve un problema de valores propios simétrico de tamaño $r \times r$.
5. Refinamiento (Opcional): Si la precisión no es suficiente debido a la sensibilidad de los subespacios (cuando las partes imaginarias están muy cerca), se aplican pasos de corrección basados en algoritmos de Jacobi sobre los subespacios problemáticos.

3. Contribuciones Clave

• Nueva Estrategia de Factorización: Se demuestra que para matrices normales reales, la descomposición puede lograrse resolviendo principalmente un problema de valores singulares (SVD) de tamaño $n/2$ y problemas de valores propios simétricos pequeños, evitando la reducción costosa a forma de Hessenberg de la matriz completa $A$.
• Análisis de Estabilidad y Sensibilidad: Se proporciona un análisis riguroso de cómo los errores de punto flotante afectan la precisión. Se demuestra que la precisión de los vectores de Schur depende de la separación (gap) entre las partes imaginarias de los valores propios.
• Complejidad Computacional:
  • En el caso ideal (matrices ortogonales aleatorias de Haar, donde casi todos los valores propios son complejos y distintos), el costo es equivalente a una factorización de Hessenberg ($\sim \frac{14}{3}n^3$ operaciones), pero sin la necesidad de iteraciones QR posteriores costosas para converger a la forma diagonal.
  • El método es significativamente más rápido que los algoritmos estándar (como gees de LAPACK) para matrices ortogonales, superándolos por un factor de 2 a 10 dependiendo del tamaño.
• Aplicación en Geometría Riemanniana: Se destaca la utilidad del método para calcular el logaritmo de matrices ortogonales y el baricentro Riemanniano (media de Karcher) en grupos de Lie como $SO(n)$, donde la velocidad es crítica y la precisión de máquina no siempre es necesaria.

4. Resultados Numéricos

Los autores validan el método mediante experimentos extensos:

• Precisión: En la mayoría de los casos (valores propios bien separados), el método alcanza una precisión cercana a la de la máquina ($\epsilon_m$). En casos de valores propios agrupados, el algoritmo de refinamiento (paso 8) restaura la precisión deseada.
• Rendimiento:
  • Para matrices ortogonales de tamaño $n=100$, el nuevo método es ~10 veces más rápido que gees (LAPACK).
  • Para $n=10000$, es ~2 veces más rápido.
  • El costo es comparable a la reducción a Hessenberg, pero evita la fase de iteración QR que suele ser el cuello de botella en los métodos tradicionales para obtener vectores propios.
• Escenarios de Peor Caso: El método puede volverse más costoso si hay grandes clústeres de partes imaginarias (requiriendo resolver EVPs grandes en el paso 5), pero esto se considera un caso poco probable en aplicaciones prácticas como matrices aleatorias.

5. Significado e Impacto

Este trabajo es significativo por varias razones:

1. Eficiencia Específica: Ofrece una solución óptima para una clase de matrices (normales reales) que es muy común en aplicaciones de física, estadística y geometría, pero que a menudo se trata con algoritmos genéricos ineficientes.
2. Explotación de Estructura: Demuestra que la estructura antisimétrica de las matrices normales es una vía más eficiente que la simétrica o la reducción general a Hessenberg, especialmente cuando el número de valores propios reales es bajo ($r \ll n$).
3. Herramienta para Manifold Learning: Proporciona un componente computacional crítico (EVD rápida y estable) para algoritmos de optimización en variedades Riemannianas, como el cálculo de medias y gradientes en el grupo ortogonal, permitiendo iteraciones más rápidas en algoritmos de descenso de gradiente.
4. Reproducibilidad: Los autores han hecho público el código en Julia (NormalEVP.jl), facilitando la adopción y verificación de sus resultados.

En resumen, el artículo presenta un algoritmo robusto y rápido que reemplaza los métodos estándar para la descomposición de matrices normales reales, aprovechando la teoría de grupos de invarianza y la descomposición de la parte antisimétrica para lograr ganancias significativas en rendimiento sin sacrificar la estabilidad numérica.