Numerically stable evaluation of closed-form expressions for eigenvalues of $3 \times 3$ matrices

Cet article propose une méthode d'évaluation fermée et numériquement stable des valeurs propres de matrices réelles diagonalisables $3 \times 3$ via quatre invariants, démontrant à la fois une précision accrue et une vitesse de calcul dix fois supérieure à la bibliothèque LAPACK pour des cas difficiles.

L'essentiel

Imagine que vous êtes un architecte ou un ingénieur qui doit construire un pont. Pour que le pont soit solide, vous devez connaître les points de pression les plus forts et les plus faibles. En mathématiques, pour les structures complexes (représentées par des grilles de nombres appelées matrices), ces points de pression sont appelés valeurs propres.

Le problème, c'est que pour les grilles de 3x3 (comme un petit carré de 9 nombres), les formules mathématiques classiques pour trouver ces points sont comme un château de cartes : très fragiles. Si deux points de pression sont presque identiques (ce qui arrive souvent dans la réalité), le château s'effondre et les calculs deviennent faux.

Voici l'explication de ce papier de recherche, racontée comme une histoire de réparation de cette fragilité.

1. Le Problème : Le "Château de Cartes" Mathématique

Dans le monde réel, les ingénieurs utilisent souvent des formules trigonométriques (avec des cosinus et des sinus) pour trouver ces valeurs propres. C'est comme utiliser une règle pour mesurer une distance.

• Le souci : Quand les valeurs sont très proches les unes des autres (comme deux notes de musique presque identiques), la "règle" classique commence à trembler. Les erreurs d'arrondi (les petits défauts inévitables des ordinateurs) s'accumulent et donnent un résultat complètement faux. C'est ce qu'on appelle l'instabilité numérique.

2. La Solution : Une Nouvelle "Boîte à Outils"

Les auteurs de ce papier, Michal Habera et Andreas Zilian, ont décidé de ne pas jeter la règle, mais de la renforcer. Ils ont créé une nouvelle méthode pour calculer ces valeurs propres de manière numériquement stable.

Pour y arriver, ils ne regardent pas les 9 nombres de la matrice individuellement. Au lieu de cela, ils utilisent 4 indicateurs magiques (qu'ils appellent des "invariants") qui résument la matrice :

1. La Trace ($I_1$) : La somme des nombres sur la diagonale. C'est comme le "poids total" de la structure.
2. Les Invariants Déviatoriques ($J_2$ et $J_3$) : Ce sont des mesures de la "forme" ou de la déformation de la structure, une fois qu'on a enlevé le poids total. Imaginez que vous enlevez la gravité pour voir comment la forme se déforme.
3. Le Discriminant ($\Delta$) : C'est le gardien de la porte. Il vous dit si les valeurs propres sont distinctes ou si elles sont en train de se coller les unes aux autres (comme des aimants qui s'attirent).

3. L'Innovation : Comment ils ont rendu les calculs solides

Le papier propose des algorithmes (des recettes de cuisine mathématiques) pour calculer ces 4 indicateurs sans se tromper, même dans les cas les plus difficiles.

• L'analogie de la balance : Imaginez que vous devez peser deux objets très lourds qui sont presque identiques. Si vous les mettez sur une balance classique, le résultat peut être nul à cause d'un petit tremblement. Les auteurs proposent une balance spéciale qui mesure la différence entre les objets avant de les peser, évitant ainsi les erreurs.
• Le discriminant : Au lieu de calculer une grosse soustraction qui annule tout (ce qui est dangereux), ils décomposent le calcul en une somme de petits produits. C'est comme construire un mur avec des briques bien ajustées plutôt que d'essayer de soulever un bloc de béton géant d'un seul coup.

4. Les Résultats : Plus rapide et plus précis

Les auteurs ont testé leur méthode contre les standards de l'industrie (des bibliothèques logicielles très puissantes comme LAPACK, utilisées par les supercalculateurs).

• Précision : Leur méthode donne des résultats exacts, même quand les valeurs propres sont presque identiques, là où les méthodes anciennes échouent.
• Vitesse : C'est là que ça devient impressionnant. Leur méthode est 10 fois plus rapide que les méthodes classiques pour ce type de calcul.
  • Pourquoi ? Les méthodes classiques utilisent des boucles qui tournent encore et encore (comme un tournevis qui tourne pour serrer une vis). La nouvelle méthode est une formule directe (comme un coup de marteau précis). Pas besoin de tourner, on tape une fois et c'est fini.

5. Pourquoi c'est important pour vous ?

Vous ne voyez peut-être pas ces calculs, mais ils sont partout :

• En génie civil : Pour savoir si un pont ou un bâtiment va s'effondrer sous le vent ou un tremblement de terre.
• En médecine : Pour analyser les images IRM du cerveau (comment l'eau circule dans les tissus).
• En mécanique : Pour simuler comment une voiture réagit à un choc.

Si les calculs sont faux à cause d'une instabilité numérique, un ingénieur pourrait penser qu'un pont est sûr alors qu'il est en danger, ou inversement.

En résumé

Ce papier est comme un manuel de réparation pour les mathématiques de l'ingénierie. Il dit : "Arrêtez d'utiliser les vieilles formules fragiles qui cassent quand les choses se ressemblent trop. Utilisez notre nouvelle boîte à outils avec 4 indicateurs robustes. Vous aurez des résultats plus précis, et vous le ferez 10 fois plus vite."

C'est une victoire pour la précision et l'efficacité, permettant aux ordinateurs de résoudre des problèmes complexes de manière plus fiable, du design de voitures à la modélisation du climat.

Résumé technique

Voici un résumé technique détaillé du papier de recherche intitulé « Numerically stable evaluation of closed-form expressions for eigenvalues of 3 × 3 matrices » (Évaluation numériquement stable d'expressions en forme fermée pour les valeurs propres de matrices 3 × 3).

1. Problématique

Le calcul des valeurs propres de matrices réelles $3 \times 3$diagonalisables repose souvent sur des formules trigonométriques classiques dérivées des solutions algébriques de l'équation cubique (travaux de Cardan et Viète). Ces formules utilisent quatre invariants : la trace ($I_1$), les invariants déviateurs ($J_2, J_3$) et le discriminant ($\Delta$).

Cependant, ces approches en forme fermée souffrent d'une instabilité numérique sévère dans l'arithmétique à précision finie, particulièrement lorsque :

• Les valeurs propres sont multiples (coalescence).
• Les valeurs propres sont très proches les unes des autres.
• La base propre est mal conditionnée (matrices non normales).

Les méthodes itératives standards (comme l'algorithme QR implémenté dans LAPACK) sont stables mais coûteuses en temps de calcul et ne permettent pas toujours de garantir un ordre spécifique des valeurs propres ou une différentiation symbolique directe. L'objectif de ce travail est de combler le fossé entre la rapidité des formules fermées et la stabilité des méthodes itératives.

2. Méthodologie

Les auteurs proposent une analyse rigoureuse de la stabilité numérique et développent de nouveaux algorithmes pour évaluer les invariants clés de manière stable.

A. Analyse de la conditionnement et des erreurs

• Utilisation du modèle IEEE 754 pour analyser les erreurs d'arrondi.
• Définition de la stabilité arrière (backward stability) et de la stabilité avant (forward stability).
• Dérivation de bornes d'erreur avant (forward error bounds) basées sur le nombre de conditionnement de la base propre ($\kappa_2(U)$) et les dérivées partielles (Jacobian) des invariants.
• Introduction d'une notion de stabilité arrière relative déviatorique, cruciale pour les invariants $J_2$ et $J_3$ qui s'annulent pour les matrices proportionnelles à l'identité.

B. Algorithmes proposés

Les auteurs rejettent les expressions polynomiales naïves (sommes de monômes) qui entraînent des annulations catastrophiques. Ils proposent des algorithmes basés sur des différences diagonales et des produits de sommes :

1. Invariant $I_1$ (Trace) : Calculé par somme directe (stable).
2. Invariant $J_2$ : Au lieu de calculer $J_2 = \frac{1}{2}\text{tr}(\text{dev}(A)^2)$, l'algorithme utilise les différences diagonales ($d_i = A_{ii} - A_{jj}$) et les produits hors-diagonaux. Cela garantit que $J_2 = 0$ exactement pour les matrices scalaires, éliminant les erreurs d'arrondi près de zéro.
3. Invariant $J_3$ (Déterminant du déviateur) : Développé en termes de différences diagonales et de produits mixtes pour éviter les erreurs lors de la singularité du déviateur.
4. Discriminant $\Delta$ : Au lieu de la formule instable $\Delta = 4J_2^3 - 27J_3^2$, l'algorithme utilise une formule somme de produits (sum-of-products) dérivée de la factorisation de Parlett. Cette méthode exprime $\Delta$ comme une somme de termes qui s'annulent simultanément lorsque la matrice approche d'une configuration à valeurs propres multiples, évitant ainsi la soustraction de deux grandes quantités proches.
5. Calcul des valeurs propres : Utilisation de la formule trigonométrique avec l'angle triple $\phi = \arctan(\dots)$ plutôt que $\arccos$, car l'arctangente est plus stable autour de zéro (cas des valeurs propres multiples) que l'arccosinus autour de 1.

3. Résultats Principaux

A. Benchmarks de Stabilité Numérique

Les auteurs ont testé leurs algorithmes sur des matrices générées avec des valeurs propres proches (chemins limites vers des valeurs propres doubles ou triples) et des transformations $U$ de conditionnement variable (orthogonal, bien conditionné, mal conditionné).

• $J_2$ : L'algorithme proposé est exactement stable et satisfait les bornes d'erreur théoriques, même pour des bases propres mal conditionnées. Il surpasse largement les méthodes "naïves" et "tensorielles" (basées sur NumPy).
• $J_3$ et $\Delta$ : Les algorithmes proposés sont stables pour les bases propres bien conditionnées ou orthogonales. Cependant, pour des bases propres très mal conditionnées, l'erreur dépasse légèrement les bornes théoriques, bien que restant bien inférieure aux méthodes naïves.
• Valeurs propres : Pour les matrices à base propre bien conditionnée, la méthode proposée atteint une précision proche de la précision machine ($\approx 10^{-16}$), comparable à LAPACK. Pour les cas mal conditionnés, LAPACK (DGEEV) reste plus précis, mais la méthode proposée reste robuste.

B. Benchmarks de Performance

• Sur un cas de test difficile (valeur propre presque double), l'algorithme proposé est environ 10 fois plus rapide que la routine LAPACK DGEEV optimisée.
• Le temps d'exécution moyen est d'environ 38 ns contre 396 ns pour LAPACK.
• La nature en forme fermée permet l'inlining du code, réduisant encore l'overhead dans les applications critiques.

C. Application : Fonction de Yield Mohr-Coulomb

L'application pratique dans la mécanique des sols (critère de plasticité Mohr-Coulomb) démontre l'impact critique de la stabilité.

• Les méthodes naïves génèrent des erreurs massives ($\approx 10^{-8}$) près des angles de Lode où les contraintes principales coalescent, rendant la prédiction de la rupture imprécise.
• L'algorithme proposé maintient une erreur absolue proche de la précision machine sur tout le spectre des états de contrainte.

4. Contributions Clés

1. Analyse théorique rigoureuse : Première dérivation de bornes d'erreur avant strictes pour les invariants $J_2, J_3$ et $\Delta$ dans le cas général de matrices non symétriques, incluant l'analyse du conditionnement de la base propre.
2. Algorithmes stables : Proposition d'algorithmes spécifiques utilisant des différences diagonales et des formules somme de produits pour éviter les annulations catastrophiques, en particulier pour $J_2$ et $\Delta$.
3. Librairie Open Source : Développement de la librairie eig3x3 (C/Python) fournissant des implémentations optimisées et stables pour les matrices symétriques et générales.
4. Preuve de concept pratique : Démonstration que les solutions en forme fermée peuvent rivaliser avec les méthodes itératives en termes de précision tout en étant nettement plus rapides, à condition que la base propre ne soit pas excessivement mal conditionnée.

5. Signification et Conclusion

Ce travail démontre que les formules fermées pour les valeurs propres de matrices $3 \times 3$ ne sont pas intrinsèquement instables, mais que leur évaluation naïve l'est. En reformulant le calcul des invariants via des techniques de stabilité numérique avancées, il est possible d'obtenir des solutions rapides, précises et différentiables.

Cela ouvre la voie à une utilisation généralisée de ces méthodes dans des domaines exigeants en performance et en dérivées symboliques, tels que :

• L'apprentissage automatique (gradients de fonctions de matrices).
• La mécanique des solides et des fluides (calcul de contraintes principales, stabilité des vortex).
• L'imagerie médicale (DTI).
• L'optimisation structurelle.

Bien que les méthodes itératives (LAPACK) restent la référence absolue pour les matrices extrêmement mal conditionnées, l'approche proposée offre une alternative robuste et ultra-rapide pour la grande majorité des cas pratiques rencontrés en ingénierie.