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

Résumé technique : La méthode d'Arnoldi par décalage et inversion pour les faisceaux de matrices singuliers

Énoncé du problème
L'article aborde le défi de la résolution de problèmes de valeurs propres généralisées singulières, grandes et creuses, de la forme $Ax = \lambda Bx$, où le faisceau de matrices $A - \lambda B$ est singulier (c'est-à-dire que $\det(A - \lambda B) \equiv 0$). Dans de tels cas, la définition classique des valeurs propres échoue car tout scalaire satisfait la condition du déterminant. Les solutions significatives sont les valeurs propres « régulières », définies là où le rang de $A - \lambda_0 B$ chute en dessous du rang normal du faisceau.

Bien que des méthodes telles que la méthode en escalier ou la décomposition QZ existent pour les problèmes petits et denses, elles sont prohibitives en termes de calcul pour les systèmes grands et creux. Des approches récentes proposées par Hochstenbach, Mehl et Plestenjak suggèrent de régulariser les faisceaux singuliers en utilisant des matrices denses choisies aléatoirement pour séparer les vraies valeurs propres des valeurs propres fantômes. Cependant, les auteurs soutiennent que l'utilisation de matrices de régularisation denses détruit la structure creuse du problème original, entraînant des surcoûts significatifs en performance et en mémoire pour les applications à grande échelle. De plus, les solveurs directs creux existants ne prennent pas nativement en charge la suppression ou l'ajout de lignes/colonnes nécessaire pour régulariser ces faisceaux sans développement personnalisé.

Méthodologie
Les auteurs proposent une stratégie déterministe pour régulariser les faisceaux singuliers en exploitant la séquence de pivotage d'une factorisation LU creuse. La méthodologie centrale comprend les étapes suivantes :

1. Régularisation par pivotage LU : Au lieu de matrices denses aléatoires, la méthode utilise une factorisation LU creuse de la matrice décalée $A - \sigma B$ (où $\sigma$ est un décalage). Au cours de la factorisation, si un pivot nul (ou proche de zéro) est rencontré, l'algorithme ajoute une ligne et une colonne à la matrice pour maintenir la factorisation.

   • Faisceau augmenté (AGEP) : Le processus construit des matrices $F$ et $G$ telles que le faisceau augmenté $\begin{bmatrix} A & G \\ F^* & 0 \end{bmatrix} - \lambda \begin{bmatrix} B & 0 \\ 0 & 0 \end{bmatrix}$ devienne régulier. Les matrices $F$ et $G$ sont dérivées directement des matrices de permutation et du processus de sélection de pivot de la factorisation LU.
   • Faisceau projeté (PGEP) : Alternativement, la méthode construit un faisceau projeté $G_\perp^* (A - \lambda B) F_\perp$ en supprimant les lignes et les colonnes correspondant aux pivots nuls. $F_\perp$ et $G_\perp$ sont dérivés de la matrice identité avec des lignes spécifiques supprimées.

2. Détection et correction du rang : La factorisation LU tente naturellement de révéler le rang normal du faisceau. Les auteurs notent que, en raison d'erreurs d'arrondi ou de paramètres de tolérance inappropriés, le rang détecté ($\hat{k}$) peut différer du vrai rang normal ($k$).

   • Si le rang détecté est trop faible (c'est-à-dire que la tolérance est trop agressive), un schéma de correction de rang est proposé. Cela implique l'utilisation de matrices orthogonales aléatoires pour augmenter les matrices $F$ et $G$ détectées à la dimension correcte, suivie d'une mise à jour peu coûteuse du complément de Schur pour résoudre le système résultant.

3. Arnoldi par décalage et inversion : Le problème régularisé (soit AGEP, soit PGEP) est résolu à l'aide de la méthode d'Arnoldi bilatérale par décalage et inversion. Les auteurs prouvent un théorème d'équivalence montrant que la factorisation d'Arnoldi sur le faisceau augmenté (avec un vecteur de départ spécifique et un produit scalaire) produit les mêmes valeurs de Ritz et la même matrice de Hessenberg que le faisceau projeté. Cela permet d'utiliser la formulation projetée pour éviter la construction explicite du grand système augmenté tout en conservant les mêmes propriétés spectrales.

4. Problèmes rectangulaires : La méthode est étendue aux problèmes de valeurs propres rectangulaires ($A, B \in \mathbb{R}^{n \times m}$ avec $n \neq m$). Pour les problèmes rectangulaires de plein rang, les solveurs creux existants avec pivotage de lignes peuvent être utilisés directement pour sélectionner un sous-problème carré et régulier.

Contributions clés

• Régularisation creuse : L'article introduit une méthode pour générer des matrices de régularisation creuses ($F, G, F_\perp, G_\perp$) directement à partir de la séquence de pivotage des factorisations LU creuses, préservant ainsi la structure creuse du problème original.
• Détection du rang : Il démontre que les factorisations LU creuses standard (utilisant un pivotage partiel, de type « rook » ou complet) détectent souvent avec succès le rang normal du faisceau, facilitant la construction automatique de la régularisation.
• Correction du rang : Un algorithme spécifique est proposé pour corriger la régularisation si la factorisation LU sous-estime le rang normal, garantissant que la méthode reste robuste même avec une détection de rang imparfaite.
• Équivalence théorique : L'article établit un lien théorique rigoureux entre les faisceaux régularisés augmentés et projetés dans le contexte de la méthode d'Arnoldi, prouvant qu'ils produisent des approximations spectrales identiques dans des conditions spécifiques.

Résultats
Des expériences numériques ont été menées sur trois types de problèmes :

1. Mise à jour de modèle par éléments finis : Un modèle de structure de treillis avec 33 degrés de liberté (aboutissant à un faisceau de dimension ~35 937). La méthode proposée (utilisant le pivotage RLU) a correctement identifié le rang et calculé les valeurs propres avec des résidus de l'ordre de $10^{-14}$. En revanche, l'approche randomisée a entraîné des nombres de conditionnement plusieurs ordres de grandeur plus élevés ($10^6$ contre $10^3$) et des résidus plus importants.
2. Détection de valeurs propres doubles : Un problème paramétrique impliquant la discrétisation d'un Laplacien 2D. La méthode basée sur LU a correctement déterminé le rang et atteint une haute précision. L'approche randomisée a nécessité la construction de matrices denses, entraînant des coûts de stockage et de calcul nettement plus élevés.
3. Problème de valeurs propres non linéaires rectangulaire : Un problème rectangulaire de plein rang dérivé d'un problème de valeurs propres non linéaire. L'approche LU creuse a été significativement plus rapide (0,05 s contre 2,9 s pour la factorisation et la projection) et a nécessité moins de mémoire (98 Mo contre 204 Mo) par rapport à l'approche dense randomisée.

Dans tous les cas, les méthodes de régularisation creuse ont préservé la structure creuse des matrices, conduisant à des nombres de conditionnement plus faibles pour les systèmes décalés et à une convergence plus rapide par rapport à la régularisation dense randomisée.

Signification et affirmations
Les auteurs positionnent ce travail comme un concept plutôt que comme une solution logicielle complète et optimisée, reconnaissant que des solveurs directs creux robustes pour les matrices singulières font actuellement défaut dans les bibliothèques logicielles standard. Cependant, l'article affirme que :

• Performance : Préserver la structure creuse via un pivotage déterministe est bénéfique tant pour la performance du solveur de valeurs propres que pour la précision des résultats (via un meilleur conditionnement).
• Faisabilité : Les solveurs directs creux existants peuvent être adaptés pour résoudre ces problèmes, en particulier pour les cas rectangulaires de plein rang, sans nécessiter de nouvelle infrastructure logicielle complexe.
• Robustesse : Le schéma de correction de rang proposé assure que la méthode fonctionne même lorsque la factorisation LU ne parvient pas à détecter parfaitement le rang normal.

L'article conclut que, bien que l'approche randomisée soit théoriquement fondée, l'approche déterministe LU creuse offre une alternative pratique, efficace et précise pour les problèmes de valeurs propres singuliers à grande échelle, à condition que l'écart entre les valeurs singulières soit suffisant pour la détection du rang. Des travaux futurs sont suggérés pour explorer des solveurs avec pivotage de volume maximal global et pour exploiter les structures de Kronecker dans les problèmes multiparamétriques.