A Jacobi-like algorithm for normal matrices by the skew-symmetric part

Résumé technique : Un algorithme de type Jacobi pour les matrices normales par la partie antisymétrique

Énoncé du problème
L'article traite du problème des valeurs propres réelles normales (EVP), qui consiste à calculer les valeurs propres et, le cas échéant, les vecteurs propres d'une matrice normale réelle $A$. Bien que l'algorithme de Jacobi soit une méthode classique pour les matrices symétriques, son extension aux matrices normales est difficile. Les généralisations existantes, telles que la méthode de Goldstine et Horwitz, reposent souvent sur l'arithmétique complexe, ce qui est plus coûteux en calcul et ignore la structure inhérente des matrices réelles. D'autres approches en arithmétique réelle, comme celles de Zhou et Brent ou de Bunse-Gerstner et al., dépendent de sous-routines pour la décomposition de Schur réelle de matrices arbitraires $4 \times 4$, ce qui peut être très intensif en calcul. Les auteurs notent que, si les algorithmes de type Jacobi ont retrouvé de la popularité grâce à leur parallélisabilité, ils souffrent souvent d'une convergence plus lente ou d'un nombre d'opérations plus élevé par rapport à des concurrents comme l'algorithme QR, en particulier pour les matrices ayant des valeurs propres complexes.

Méthodologie
La méthode proposée est un algorithme de type Jacobi qui exploite la commutativité des parties symétrique et antisymétrique d'une matrice normale réelle. L'innovation centrale réside dans l'utilisation de la méthode de Paardekooper — un algorithme rapide à convergence quadratique pour les matrices antisymétriques — comme moteur principal de la décomposition.

L'algorithme (Algorithme 5.1) se déroule en trois étapes distinctes :

1. Étape I : Méthode implicite de Paardekooper.
   L'algorithme applique d'abord implicitement la méthode de Paardekooper à la partie antisymétrique $\Omega = \text{skew}(A)$. Au lieu de transformer $\Omega$ directement, les rotations dérivées de $\Omega$ sont appliquées à la matrice complète $A$. Cette étape se poursuit jusqu'à ce que les blocs hors diagonale de $\Omega$ soient réduits en dessous d'un seuil prescrit $\rho$. Théoriquement, cela rend la partie antisymétrique bloc-diagonale, ce qui, pour les matrices normales, révèle la structure des valeurs propres.

2. Étape II : Identification des blocs et diagonalisation spécialisée.
   Après l'Étape I, la matrice $A$ n'est pas encore entièrement bloc-diagonale. L'algorithme identifie les composantes connexes (sous-matrices non nulles) dans $A$ basées sur le seuil $\sqrt{\rho}\|A\|_F$. Ces composantes sont catégorisées en trois types selon une analyse de perturbation :

   • Blocs Symétriques Skew-Hamiltoniens (SSkH) : Si un bloc correspond à des valeurs propres avec des parties imaginaires non nulles répétées, il est traité comme une matrice SSkH perturbée. L'algorithme calcule la matrice SSkH la plus proche et la diagonalise en utilisant une méthode spécialisée de type Jacobi (Algorithme 5.3) adaptée des techniques pour les matrices symétriques skew-Hamiltoniennes.
   • Blocs Symétriques : Si un bloc correspond à des valeurs propres réelles, il est traité comme une matrice symétrique perturbée. L'algorithme applique la méthode de Jacobi standard pour les matrices symétriques (Algorithme 5.4) à la partie symétrique du bloc.
   • Blocs Généraux : Si la perturbation est trop importante pour les méthodes spécialisées (par exemple, des valeurs propres avec des parties imaginaires proches mais distinctes), le bloc est transmis à un solveur général pour matrices normales.

3. Étape III : Raffinement de la précision.
   Une étape finale applique un algorithme standard de type Jacobi pour les matrices normales (spécifiquement, la méthode de Zhou et Brent) à l'ensemble de la matrice. Parce que les Étapes I et II ont déjà réduit considérablement les éléments hors diagonale, cette étape finale bénéficie d'une convergence locale quadratique, nécessitant très peu de balayages (généralement un ou deux) pour atteindre la précision cible $\rho$.

Contributions clés

• Conception algorithmique : L'article introduit un algorithme hybride de type Jacobi qui exploite la rapidité de la méthode de Paardekooper pour les matrices antisymétriques afin d'accélérer la résolution de l'EVP réel normal.
• Transformation implicite : La méthode applique les transformations dérivées de la partie antisymétrique directement à la matrice normale, préservant l'arithmétique réelle et évitant les opérations sur les nombres complexes.
• Analyse de perturbation : Les auteurs fournissent une analyse de perturbation rigoureuse (Section 6) qui quantifie comment l'imprécision dans la décomposition de Schur de la partie antisymétrique affecte la structure de la matrice normale. Cette analyse justifie la logique de branchement spécifique de l'Étape II, déterminant quand utiliser des solveurs spécialisés SSkH ou symétriques par rapport à un solveur général.
• Formules de matrices les plus proches : L'article dérive des formules explicites pour la matrice skew-Hamiltonienne symétrique la plus proche et la matrice ortho-symplectique la plus proche (Section 8), qui sont cruciales pour les bornes de perturbation et la conception des étapes de diagonalisation spécialisées.

Résultats
Les expériences numériques (Section 7) comparent l'Algorithme 5.1 à quatre autres algorithmes de type Jacobi : Bunse-Gerstner et al., Goldstine et Horwitz, Zhou et Brent, et une méthode randomisée (RandDiag).

• Précision : L'Algorithme 5.1 atteint systématiquement la précision cible ($\text{offschur}(A) \leq \rho \|A\|_F$) sur diverses distributions de matrices normales, y compris celles avec des valeurs propres réelles, des parties imaginaires répétées et des matrices orthogonales aléatoires. Il se comporte de manière comparable aux concurrents les plus précis.
• Efficacité : La méthode proposée démontre un avantage significatif en vitesse, étant 5 à 10 fois plus rapide que les autres algorithmes de type Jacobi en termes de temps d'exécution. Ce gain est attribué à la convergence rapide de la méthode de Paardekooper dans les étapes initiales, ce qui réduit la charge de travail pour l'étape de raffinement finale.
• Robustesse : La méthode reste robuste même lorsque les valeurs propres ont des parties imaginaires répétées ou très proches, des scénarios où d'autres méthodes pourraient avoir des difficultés ou nécessiter plus d'itérations.

Signification
L'article affirme que l'algorithme proposé offre une alternative computationnellement supérieure pour résoudre le problème des valeurs propres réelles normales, en particulier pour les matrices où la plupart des valeurs propres sont complexes, telles que les matrices orthogonales aléatoires apparaissant en statistique sur les variétés. En combinant les propriétés structurelles des matrices normales avec l'efficacité de la méthode de Paardekooper, l'algorithme atteint une haute parallélisabilité et rapidité sans sacrifier la précision. Les auteurs positionnent ce travail comme une amélioration pratique par rapport aux variantes existantes de type Jacobi, notant que la méthode est la plus efficace dans les contextes où la partie antisymétrique fournit un guide structurel fort pour la décomposition. L'article conclut en soulignant la disponibilité d'une implémentation open-source pour la reproductibilité.