Méthode de la puissance inverse

En analyse numérique, la méthode de puissance inverse est un algorithme de recherche de valeur propre itératif de valeurs propres. Il permet de trouver un vecteur propre approché lorsqu'une approximation d'une valeur propre correspondante est déjà connue. La méthode est conceptuellement similaire à la méthode de la puissance itérée. Il semble avoir été initialement développé pour calculer les fréquences de résonance dans le domaine de la mécanique des structures.

La méthode de la puissance inverse commence par une approximation $\mu$ pour la valeur propre correspondant au vecteur propre souhaité et un vecteur $b_{0}$, soit un vecteur sélectionné au hasard, soit une approximation du vecteur propre. La méthode est décrite par la récurrence :

$$b_{k+1}={\frac {(A-\mu I)^{-1}b_{k}}{C_{k}}},$$

où les $C_{k}$ sont des constantes habituellement choisies telles que $C_{k}=\|(A-\mu I)^{-1}b_{k}\|.$ Puisque les vecteurs propres sont définis à une constante multiplicative près, le choix de $C_{k}$ peut être arbitraire en théorie ; les aspects pratiques du choix de $C_{k}$ sont discutés ci-dessous.

À chaque itération, le vecteur $b_{k}$ est multiplié par la matrice $(A-\mu I)^{-1}$ et normalisé. On retrouve la même formule que dans la méthode de la puissance itérée, mais en remplaçant la matrice $A$ par $(A-\mu I)^{-1}.$ Plus la valeur choisie $\mu$ est proche de la valeur propre réelle, plus l'algorithme converge rapidement ; cependant, un choix incorrect de $\mu$ peut conduire à une convergence lente ou à la convergence vers un vecteur propre différent de celui souhaité. En pratique, la méthode est utilisée lorsqu’une bonne approximation de la valeur propre est connue et qu’il suffit donc de quelques itérations (souvent une seule).

Théorie et convergence

L'idée de base de la puissance itérée est de choisir un vecteur initial $b$ (soit une approximation de vecteur propre, soit un vecteur aléatoire) et de calculer itérativement $Ab,A^{2}b,A^{3}b,...$. Sauf pour un ensemble de mesure nulle, pour tout vecteur initial, le résultat convergera vers un vecteur propre correspondant à la valeur propre dominante de la matrice.

La puissance inverse fait de même pour la matrice $(A-\mu I)^{-1}$, donc il converge vers le vecteur propre correspondant à la valeur propre dominante de la matrice $(A-\mu I)^{-1}$. Les valeurs propres de cette matrice sont $(\lambda _{1}-\mu )^{-1},...,(\lambda _{n}-\mu )^{-1},$ où $\lambda _{i}$ sont les valeurs propres de $A$. Le plus grand de ces nombres correspond au plus petit nombre parmi $(\lambda _{1}-\mu ),...,(\lambda _{n}-\mu ).$ Les vecteurs propres de $A$ et de $(A-\mu I)^{-1}$ sont les mêmes, puisque

$$Av=\lambda v\Leftrightarrow (A-\mu I)v=\lambda v-\mu v\Leftrightarrow (\lambda -\mu )^{-1}v=(A-\mu I)^{-1}v$$

Ainsi, la méthode converge vers le vecteur propre de la matrice $A$ correspondant à la valeur propre la plus proche de $\mu .$

En particulier, en prenant $\mu =0$ on voit que $A^{-1}b_{k}$ converge vers le vecteur propre correspondant à la valeur propre de $A^{-1}$ avec la plus grande valeur absolue ${\frac {1}{\lambda _{N}}}$ et peut donc être utilisé pour déterminer la plus petite valeur propre de grandeur de $A$ puisqu'ils sont inversement liés.

Vitesse de convergence

On détermine ici le taux de convergence de la méthode.

On sait que l'algorithme de la puissance itérée converge linéairement vers la limite, plus précisément :

$$\mathrm {d} (b^{\mathrm {ideal} },b_{\mathrm {puiss.it.} }^{k})=O\left(\left|{\frac {\lambda _{\min }}{\lambda _{\max }}}\right|^{k}\right),$$

par conséquent, pour la méthode de la puissance inverse, un résultat similaire ressemble à :

$$\mathrm {d} (b^{\mathrm {ideal} },b_{\mathrm {puiss.inv.} }^{k})=O\left(\left|{\frac {\mu -\lambda _{\mathrm {plus\ proche\ de\ } \mu }}{\mu -\lambda _{2^{e}\ \mathrm {plus\ proche\ de\ } \mu }}}\right|^{k}\right).$$

Il s'agit d'une formule clé pour comprendre la convergence de la méthode. Cela montre que si $\mu$ est choisi suffisamment proche d'une valeur propre $\lambda$, Par exemple $\mu -\lambda =\epsilon$ chaque itération améliorera la précision $|\epsilon |/|\lambda +\epsilon -\lambda _{{\text{plus proche de}}\lambda }|$ fois. (on l'utilise pour des choses assez petites $\epsilon$ "le plus proche de $\mu$ " et " le plus proche de $\lambda$ " est la même chose.) Pour assez petit $|\epsilon |$ c'est à peu près la même chose que $|\epsilon |/|\lambda -\lambda _{{\text{plus proche de}}\lambda }|$. Donc si l'on parvient à trouver $\mu$, de telle sorte que le $\epsilon$ sera suffisamment petit, alors très peu d'itérations peuvent être satisfaisantes.

Complexité

L'algorithme d'itération inverse nécessite la résolution d'un système linéaire ou le calcul de la matrice inverse. Pour les matrices non structurées (ni creuses, ni de Toeplitz,...) cela nécessite $O(n^{3})$ opérations.

Options de mise en œuvre

La méthode est définie par la formule :

$$b_{k+1}={\frac {(A-\mu I)^{-1}b_{k}}{C_{k}}},$$

Il existe cependant plusieurs options pour sa mise en œuvre.

Calculer une matrice inverse ou résoudre un système d'équations linéaires

On peut réécrire la formule de la manière suivante :

$$(A-\mu I)b_{k+1}={\frac {b_{k}}{C_{k}}},$$

soulignant que pour trouver l'approximation suivante $b_{k+1}$ on peut résoudre un système d'équations linéaires. Il existe deux options : on peut choisir un algorithme qui résout un système linéaire, ou on peut calculer l'inverse $(A-\mu I)^{-1}$ puis l'appliquer au vecteur. Les deux options ont une complexité O( n³), le nombre exact dépend de la méthode choisie.

Le choix dépend également du nombre d'itérations. Naïvement, si à chaque itération on résout un système linéaire, la complexité sera k O( n³), où k est le nombre d'itérations ; de même, calculer la matrice inverse et l'appliquer à chaque itération est de complexité k O( n³). Notons cependant que si l’estimation des valeurs propres $\mu$ reste constant, alors on peut réduire la complexité à O( n³) + k O(n ²) avec l'une ou l'autre méthode. Calculer la matrice inverse une fois et la stocker pour l'appliquer à chaque itération est de complexité O( n³) + k O(n ²). Stocker une décomposition LU de $(A-\mu I)$ et utiliser l'algorithme de montée et descente pour résoudre le système d'équations à chaque itération est également de complexité O( n³) + k O(n ²).

L'inversion de la matrice aura généralement un coût initial plus élevé, mais un coût inférieur à chaque itération. À l’inverse, la résolution de systèmes d’équations linéaires aura généralement un coût initial moindre, mais nécessitera plus d’opérations pour chaque itération.

Tridiagonalisation, forme de Hessenberg

S'il est nécessaire d'effectuer de nombreuses itérations (ou peu d'itérations, mais pour de nombreux vecteurs propres), alors il pourrait être judicieux de mettre d'abord la matrice sous forme de Hessenberg supérieure (pour une matrice symétrique, ce sera une forme tridiagonale), qui a un coût en ${\textstyle {\frac {10}{3}}n^{3}+O(n^{2})}$ opérations arithmétiques utilisant une technique basée sur la transformation de Householder, avec une suite finie de transformations de similarité orthogonales, un peu comme une décomposition QR bilatérale. (Pour la décomposition QR, les rotations de Householder sont multipliées uniquement à gauche, mais pour le cas de Hessenberg, elles sont multipliées à gauche et à droite). Pour les matrices symétriques, cette procédure coûte ${\textstyle {\frac {4}{3}}n^{3}+O(n^{2})}$ opérations arithmétiques utilisant une technique basée sur la réduction de Householder.

La solution du système d'équations linéaires pour les matrices tridiagonales est en $O(n)$ opérations, donc la complexité augmente en $O(n^{3})+kO(n)$, où $k$ est le numéro d'itération, ce qui est meilleur que pour l'inversion directe. Cependant, pour quelques itérations, une telle transformation peut ne pas être pratique.

La transformation vers une forme de Hessenberg implique également des racines carrées et l'opération de division, qui ne sont pas universellement prises en charge par le matériel.

Choix de la constante de normalisation C_k

Sur les processeurs à usage général (produits par exemple par Intel), le temps d'exécution de l'addition, de la multiplication et de la division est approximativement égal. Mais sur le matériel embarqué et/ou à faible consommation d'énergie (processeurs de signaux numériques, FPGA, ASIC), la division peut ne pas être prise en charge par le matériel et doit donc être évitée. Choisir $C_{k}=2^{n_{k}}$ permet une division rapide sans prise en charge matérielle explicite, car la division par une puissance de 2 peut être implémentée soit sous la forme d'un décalage de bits (pour l'arithmétique à virgule fixe), soit sous la forme d'une soustraction de $k$ de l'exposant (pour l'arithmétique à virgule flottante).

Lors de la mise en œuvre de l'algorithme utilisant l'arithmétique à virgule fixe, le choix de la constante $C_{k}$ est particulièrement important. De petites valeurs entraîneront une croissance rapide de la norme de $b_{k}$ et dépasser ; de grandes valeurs de $C_{k}$ va faire tendre le vecteur $b_{k}$ vers zéro.

Utilisation

La principale application de la méthode est la situation dans laquelle une approximation d’une valeur propre est trouvée et où il faut trouver le vecteur propre approché correspondant. Dans une telle situation, l’itération inverse est la méthode principale et probablement la seule à utiliser.

Méthodes pour trouver des valeurs propres approximatives

Généralement, la méthode est utilisée en combinaison avec une autre méthode qui trouve des valeurs propres approchées : l'exemple standard est l'algorithme de bissection pour les valeurs propres, un autre exemple est l'itération du quotient de Rayleigh, qui est en fait la même itération inverse avec le choix de la valeur propre approchée comme valeur propre, le quotient de Rayleigh correspondant au vecteur obtenu à l'étape précédente de l'itération.

Il existe certaines situations où la méthode peut être utilisée seule, mais elles sont assez marginales.

Norme de matrice comme approximation de la valeur propre dominante

La valeur propre dominante peut être facilement estimée pour n’importe quelle matrice. Pour toute norme matricielle induite, il est vrai que $\left\|A\right\|\geq |\lambda |,$ pour toute valeur propre $\lambda$. Ainsi, en prenant la norme de la matrice comme valeur propre approximative, on peut voir que la méthode convergera vers le vecteur propre dominant.

Estimations basées sur des statistiques

Dans certaines applications en temps réel, il faut trouver des vecteurs propres pour des matrices à une vitesse de plusieurs millions de matrices par seconde. Dans de telles applications, les statistiques des matrices sont généralement connues à l’avance et on peut prendre comme valeur propre approximative la valeur propre moyenne d’un grand échantillon matriciel. On peut aller plus loin et calculer le rapport moyen des valeurs propres à la trace ou à la norme de la matrice et estimer la valeur propre moyenne comme la trace ou la norme multipliée par la valeur moyenne de ce rapport. Il est clair qu’une telle méthode ne peut être utilisée qu’avec parcimonie et uniquement lorsqu’une grande précision n’est pas critique. Cette approche d'estimation d'une valeur propre moyenne peut être combinée avec d'autres méthodes pour éviter une erreur trop importante.

Voir également

• Méthode de la puissance itérée
• Itération du quotient de Rayleigh
• Algorithme de recherche de valeur propre

Références

• (en) Cet article est partiellement ou en totalité issu de l’article de Wikipédia en anglais intitulé « Inverse iteration » (voir la liste des auteurs).