Algorithme QR

En algèbre linéaire numérique, l' algorithme QR ou itération QR est un algorithme de recherche de valeurs propres, une procédure permettant de calculer les valeurs et vecteurs propres d'une matrice . L'algorithme QR a été développé à la fin des années 1950 par John G.F. Francis et par Vera N. Kublanovskaïa, travaillant indépendamment. L'idée de base est d'effectuer une décomposition QR, en écrivant la matrice comme un produit d'une matrice orthogonale et d'une matrice triangulaire supérieure, en multipliant les facteurs dans l'ordre inverse et en itérant.

Algorithme QR en pratique

Formellement, soit A une matrice réelle dont on veut calculer les valeurs propres, et soit A₀ := A . À la k -ième étape (en commençant par k = 0 ), on calcule la décomposition QR A_k = Q_k R_k où Q_k est une matrice orthogonale (c'est-à-dire Q^T = Q⁻¹ ) et R_k est une matrice triangulaire supérieure. On forme alors A_k+1 = R_k Q_k . On peut alors noter que $$A_{k+1}=R_{k}Q_{k}=Q_{k}^{-1}Q_{k}R_{k}Q_{k}=Q_{k}^{-1}A_{k}Q_{k}=Q_{k}^{\mathsf {T}}A_{k}Q_{k},$$ ainsi tous les A_k sont semblables et ont donc les mêmes valeurs propres. L'algorithme est numériquement stable car il procède par transformations de similarité orthogonales .

Sous certaines conditions, les matrices A _k convergent vers une matrice triangulaire, la forme de Schur de A . Les valeurs propres d'une matrice triangulaire sont ses éléments diagonaux et le problème des valeurs propres est résolu. Lors des tests de convergence, il n’est pas pratique d’exiger des zéros exacts,^{[réf. nécessaire]} mais le théorème de Gershgorin fournit une limite à l'erreur.

En utilisant la forme de Hessenberg

Dans la forme brute ci-dessus, les itérations sont relativement coûteuses. Cela peut être atténué en prenant d'abord la forme de Hessenberg supérieure de la matrice A (ce qui coûte ${\textstyle {\tfrac {10}{3}}n^{3}+{\mathcal {O}}(n^{2})}$ opérations arithmétiques utilisant une technique basée sur la réduction 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 réflecteurs de Householder sont multipliés uniquement à gauche, mais pour le cas de Hessenberg, ils sont multipliés à la fois à gauche et à droite.) Déterminer la décomposition QR d'une matrice de Hessenberg supérieure coûte ${\textstyle 6n^{2}+{\mathcal {O}}(n)}$ opérations arithmétiques. De plus, comme la forme de Hessenberg est déjà presque triangulaire supérieure (elle n'a qu'une seule entrée non nulle sous chaque diagonale), son utilisation comme point de départ réduit le nombre d'étapes nécessaires à la convergence de l'algorithme QR.

Si la matrice d'origine est symétrique, alors la matrice de Hessenberg supérieure est également symétrique et donc tridiagonale, et il en va de même pour toutes les A_k . Dans ce cas, atteindre la forme Hessenberg coûte ${\textstyle {\tfrac {4}{3}}n^{3}+{\mathcal {O}}(n^{2})}$ opérations arithmétiques utilisant une technique basée sur la réduction de Householder. Déterminer la décomposition QR d'une matrice tridiagonale symétrique coûte ${\mathcal {O}}(n)$ opérations.

Phase d'itération

Si une matrice de Hessenberg $A$ a un élément $a_{k,k-1}=0$ pour certains $k$, c'est-à-dire que si l'un des éléments de la sous-diagonale est nul, alors il se décompose en blocs dont les problèmes propres peuvent être résolus séparément ; une valeur propre est soit une valeur propre de la sous-matrice des premières $k-1$ lignes et colonnes, ou une valeur propre de la sous-matrice des lignes et colonnes restantes. Le but de l'étape d'itération QR est de réduire l'un de ces éléments $a_{k,k-1}$ de sorte qu'un petit bloc le long de la diagonale soit effectivement séparé de bloc principal de la matrice. Dans le cas d'une valeur propre réelle qui est généralement un bloc $1\times 1$ dans le coin inférieur droit (auquel cas l'élément $a_{nn}$ contient cette valeur propre), alors que dans le cas d'une paire de valeurs propres complexes conjuguées, c'est un bloc $2\times 2$ dans le coin inférieur droit.

Le taux de convergence dépend de la séparation entre les valeurs propres, donc un algorithme pratique utilisera des décalages, explicites ou implicites, pour augmenter la séparation et accélérer la convergence. Un algorithme QR symétrique typique isole chaque valeur propre (puis réduit la taille de la matrice) avec seulement une ou deux itérations, ce qui le rend efficace et robuste.^[pas clair]

Une seule itération avec décalage explicite

Les étapes d'une itération QR avec décalage explicite sur une matrice de Hessenberg réelle $A$ sont:

1. Choisir un décalage $\mu$ et le soustraire de tous les éléments diagonaux, produisant la matricex $A-\mu I$. Une stratégie classique est de prendre $\mu =a_{n,n}$, mais ce choix n'est pas optimal pour accélérer la convergence. Le mieux est de choisir $\mu$ proche d'une valeur propre, afin d'accélérer la convergence vers cette valeur propre.
2. Calculer la suite des rotations de Givens $G_{1},G_{2},\dots ,G_{n-1}$ sur $A-\mu I$, où$G_{i}$ agit sur les lignes $i$ et $i+1$, et $G_{i}$ est choisi pour annuler le coefficient en $(i+1,i)$ de $G_{i-1}\dotsb G_{1}(A-\mu I)$. Cela produit une matrice triangulaire supérieure $R=G_{n-1}\dotsb G_{1}(A-\mu I)$. Le facteur orthogonal $Q$ sera alors $G_{1}^{\mathrm {T} }G_{2}^{\mathrm {T} }\dotsb G_{n-1}^{\mathrm {T} }$, mais il n'est ni nécessaire ni efficace de calculer cette matrice explicitement.
3. On multiplie $R$ par les matrices de Givens $G_{1}^{\mathrm {T} }$, $G_{2}^{\mathrm {T} }$, …, $G_{n-1}^{\mathrm {T} }$ à droite, où $G_{i}^{\mathrm {T} }$ agit sur les colonnes $i$ et $i+1$. On obtient la matrice $RQ=RG_{1}^{\mathrm {T} }G_{2}^{\mathrm {T} }\dotsb G_{n-1}^{\mathrm {T} }$, qui est également une matrice de Hessenberg.
4. On annule le décalage en ajoutant $\mu$ à toutes les valeurs diagonales. Le résultat est alors $A'=RQ+\mu I$. Puisque $Q$ commute avec $I$, on retrouve $A'=Q^{\mathrm {T} }(A-\mu I)Q+\mu I=Q^{\mathrm {T} }AQ$.

Le but du décalage est de modifier les rotations Givens choisies.

Plus en détail, la structure de l'une des matrices $G_{i}$ est $$G_{i}={\begin{bmatrix}I&0&0&0\\0&c&-s&0\\0&s&c&0\\0&0&0&I\end{bmatrix}}$$ où le $I$ dans le coin supérieur gauche se trouve un bloc identité $(n-1)\times (n-1)$ et les deux scalaires $c=\cos \theta$ et $s=\sin \theta$ sont déterminés par l'angle de rotation $\theta$ approprié pour annuler le coefficient $(i+1,i)$ . Il n'est pas nécessaire d'expliciter $\theta$ ; les facteurs $c$ et $s$ peuvent être déterminés directement à partir des éléments sur lesquels la matrice $G_{i}$ devrait agir. Il n'est pas non plus nécessaire de produire la matrice entière ; multiplier (à gauche) par $G_{i}$ n'affecte que les lignes $i$ et $i+1$, il est donc plus simple de simplement mettre à jour ces deux lignes directement. De même, pour la multiplication par $G_{i}^{\mathrm {T} }$ à droite à l'étape 3, il suffit de calculer $i$, $c$, et $s$ .

Si on utilise le facteur $\mu =a_{n,n}$, alors au début de l'étape 2, on a $$A-a_{n,n}I={\begin{pmatrix}\times &\times &\times &\times &\times \\\times &\times &\times &\times &\times \\0&\times &\times &\times &\times \\0&0&\times &\times &\times \\0&0&0&\times &0\end{pmatrix}}$$ où le signe $\times$ désigne une valeur quelconque. La première rotation de Givens $G_{1}$ annuler le coefficient en $(i+1,i)$, on a alors $$G_{1}(A-a_{n,n}I)={\begin{pmatrix}\times &\times &\times &\times &\times \\0&\times &\times &\times &\times \\0&\times &\times &\times &\times \\0&0&\times &\times &\times \\0&0&0&\times &0\end{pmatrix}}{\text{.}}$$ Chaque nouvelle rotation annule un autre élément sous-diagonale, ce qui augmente le nombre de zéros connus jusqu'à ce que $$H=G_{n-2}\dotsb G_{1}(A-a_{n,n}I)={\begin{pmatrix}\times &\times &\times &\times &\times \\0&\times &\times &\times &\times \\0&0&\times &\times &\times \\0&0&0&h_{n-1,n-1}&h_{n-1,n}\\0&0&0&h_{n,n-1}&0\end{pmatrix}}{\text{.}}$$ La dernière rotation $G_{n-1}$ a des valeurs $(c,s)$ choisis telles que $sh_{n-1,n-1}+ch_{n,n-1}=0$. Si $|h_{n-1,n-1}|\gg |h_{n,n-1}|$, ce qui arrive couramment quand on arrive à la fin du processus de convergence, alors $c\approx 1$ et $|s|\ll 1$. Ainsi, la rotation donne $$R=G_{n-1}G_{n-2}\dotsb G_{1}(A-a_{n,n}I)={\begin{pmatrix}\times &\times &\times &\times &\times \\0&\times &\times &\times &\times \\0&0&\times &\times &\times \\0&0&0&\times &ch_{n-1,n}\\0&0&0&0&sh_{n-1,n}\end{pmatrix}}{\text{,}}$$ qui est la matrice triangulaire supérieure attendue. On en est alors à l'étape 3 et on effectue les rotations sur les données des colonnes. La première rotation agit sur les colonnes $1$ et $2$, produisant $$RG_{1}^{\mathrm {T} }={\begin{pmatrix}\times &\times &\times &\times &\times \\\times &\times &\times &\times &\times \\0&0&\times &\times &\times \\0&0&0&\times &ch_{n-1,n}\\0&0&0&0&sh_{n-1,n}\end{pmatrix}}{\text{.}}$$ Le motif attendu est que chaque rotation déplace une valeur non nulle de la diagonale vers la sous-diagonale, ramenant la matrice vers une forme de Hessenberg. Ceci s'arrête à $$RG_{1}^{\mathrm {T} }\dotsb G_{n-1}^{\mathrm {T} }={\begin{pmatrix}\times &\times &\times &\times &\times \\\times &\times &\times &\times &\times \\0&\times &\times &\times &\times \\0&0&\times &\times &\times \\0&0&0&-s^{2}h_{n-1,n}&csh_{n-1,n}\end{pmatrix}}{\text{.}}$$ Algébriquement la forme est inchangée, mais numériquement l'élément en position $(n,n-1)$ est devenu bien plus petit : il y a un facteur $s$ entre lui et l'élément diagonal au-dessus, mais maintenant, l'écart est plus d'un facteur $s^{2}$, et une nouvelle itération amènerait à un facteur $s^{4}$: il y a convergence quadratique. En pratique, cela signifie $O(1)$ itérations par valeur propre suffit pour avoir la convergence, et donc en tout, on peut compléter l'algorithme en $O(n)$ étapes QR, chacune d'entre elle exécutant $O(n^{2})$ opérations arithmétiques (ou au moins $O(n)$ opérations, dans le cas où $A$ est symétrique).

Visualisation

Figure 1 : Comment la sortie d'une seule itération de l'algorithme QR ou LR varie en fonction de son entrée

L'algorithme QR de base peut être visualisé dans le cas où A est une matrice symétrique définie positive. Dans ce cas, A peut être représenté comme une ellipse en 2 dimensions ou un ellipsoïde en dimensions supérieures. La relation entre l’entrée de l’algorithme et une itération unique peut alors être représentée comme dans la figure 1. Le résultat de l'algorithme LR est représenté à côté de l'algorithme QR.

Une seule itération provoque l'inclinaison ou la « chute » de l'ellipse vers l'axe des x. Dans le cas où le grand demi-axe de l'ellipse est parallèle à l'axe des x, une itération de QR ne fait rien. Une autre situation dans laquelle l'algorithme « ne fait rien » est lorsque le grand demi-axe est parallèle à l'axe des y au lieu de l'axe des x. Dans ce cas, l’ellipse peut être considérée comme étant en équilibre instable sans pouvoir tomber dans un sens ou dans l’autre. Dans les deux cas, la matrice est diagonale. Une situation dans laquelle une itération de l'algorithme « ne fait rien » est appelée un point fixe . La stratégie employée par l'algorithme est une méthode du point fixe . On observe alors qu’un point fixe est stable tandis que l’autre est instable. Si l'ellipse était inclinée très légèrement par rapport au point fixe instable, une itération de l'algorithme QR entraînerait l'inclinaison de l'ellipse à l'opposé du point fixe au lieu de s'en rapprocher. Finalement, l’algorithme convergerait vers un point fixe différent, mais cela prendrait beaucoup de temps.

Recherche de valeurs propres ou des vecteurs propres

Figure 2 : Comment la sortie d'une seule itération de QR ou LR est affectée lorsque deux valeurs propres se rapprochent l'une de l'autre

Il convient de souligner que trouver même un seul vecteur propre d'une matrice symétrique n'est pas calculable (en arithmétique réelle exacte selon les définitions de l'analyse calculable ). Cette difficulté existe chaque fois que les multiplicités des valeurs propres d'une matrice ne sont pas connaissables. En revanche, le même problème ne se pose pas pour trouver des valeurs propres. Les valeurs propres d'une matrice sont toujours calculables.

On peut discuter de la manière dont ces difficultés se manifestent dans l’algorithme QR de base. Ceci est illustré dans la figure 2. Il faut rappeler que les ellipses représentent des matrices symétriques définies positives. Lorsque les deux valeurs propres de la matrice d'entrée se rapprochent l'une de l'autre, l'ellipse d'entrée se transforme en cercle, soit un multiple de la matrice identité (un quasi-cercle correspond à un quasi-multiple de la matrice identité dont les valeurs propres sont presque égales aux entrées diagonales de la matrice). Par conséquent, le problème de la recherche approximative des valeurs propres s'avère facile dans ce cas. Cependant, une itération de QR (ou LR) incline de moins en moins les demi-axes à mesure que l'ellipse d'entrée se rapproche d'un cercle. Les vecteurs propres ne peuvent être connus que lorsque les demi-axes sont parallèles à l'axe des x et à l'axe des y. Le nombre d'itérations nécessaires pour obtenir un quasi-parallélisme augmente sans limite à mesure que l'ellipse d'entrée devient plus circulaire.

Bien qu'il puisse être impossible de calculer la décomposition propre d'une matrice symétrique arbitraire, il est toujours possible de perturber la matrice d'une quantité arbitrairement petite et de calculer la décomposition propre de la matrice résultante. Dans le cas où la matrice est représentée comme une forme quasi-circulaire, la matrice peut être remplacée par une matrice dont la représentation est un cercle parfait. Dans ce cas, la matrice est un multiple de la matrice identité et sa décomposition en valeurs propres est immédiate. Il faut alors garder à l'esprit que la base propre résultante peut être assez éloignée de la base propre d'origine.

Accélération : Déplacement et déflation

Le ralentissement lorsque l'ellipse devient plus circulaire a un effet inverse : il s'avère que lorsque l'ellipse devient plus étirée - et moins circulaire - alors la rotation de l'ellipse devient plus rapide. Un tel étirement peut être induit lorsque la matrice $M$ que l'ellipse représente est remplacé par $M-\lambda I$ où $\lambda$ est approximativement la plus petite valeur propre de $M$ . Dans ce cas, le rapport des deux demi-axes de l'ellipse se rapproche de l'infini . Dans les dimensions supérieures, un tel décalage rend la longueur du plus petit demi-axe d'un ellipsoïde petite par rapport aux autres demi-axes, ce qui accélère la convergence vers la plus petite valeur propre, mais n'accélère pas la convergence vers les autres valeurs propres. Cela devient inutile lorsque la plus petite valeur propre est entièrement déterminée, la matrice doit alors être réduite, ce qui signifie simplement supprimer sa dernière ligne et sa dernière colonne.

Le problème du point fixe instable doit également être résolu. L’heuristique de décalage est souvent conçue pour traiter également ce problème : les décalages pratiques sont souvent discontinus et aléatoires. Le décalage de Wilkinson, qui convient bien aux matrices symétriques comme celles visualisées supra, est particulièrement discontinu.

Algorithme QR implicite

Dans la pratique informatique moderne, l'algorithme QR est exécuté dans une version implicite qui facilite l'introduction de plusieurs décalages. La matrice est d'abord amenée à la forme de Hessenberg supérieure $A_{0}=QAQ^{\mathsf {T}}$ comme dans la version explicite ; puis, à chaque étape, la première colonne de $A_{k}$ est transformée via une transformation de Householder de petite taille vers la première colonne de $p(A_{k})$^[pas clair] (ou $p(A_{k})e_{1}$ ), où $p(A_{k})$, de degré $r$, est le polynôme qui définit la stratégie de changement de vitesse (souvent $p(x)=(x-\lambda )(x-{\bar {\lambda }})$, où $\lambda$ et ${\bar {\lambda }}$ sont les deux valeurs propres de la queue $2\times 2$ sous-matrice principale de $A_{k}$, le soi-disant double décalage implicite). Puis des transformations successives de taille de Householder $r+1$ sont effectuées afin de renvoyer la matrice de travail $A_{k}$ à la forme supérieure de Hessenberg. Cette opération est connue sous le nom de « bulge chasing », en raison de la forme particulière des entrées non nulles de la matrice le long des étapes de l'algorithme. Comme dans la première version, la réduction est effectuée dès qu'une des entrées sous-diagonales de $A_{k}$ est suffisamment petit.

Proposition de changement de nom

Étant donné que dans la version implicite moderne de la procédure, aucune décomposition QR n'est explicitement effectuée, certains auteurs, par exemple Watkins, ont suggéré de changer son nom en algorithme de Francis . Golub et Van Loan utilisent le terme Francis QR step .

Interprétation et convergence

L'algorithme QR peut être considéré comme une variante plus sophistiquée de l' algorithme de la puissance itérée , qui consiste à multiplier à plusieurs reprises A par un seul vecteur, en normalisant après chaque itération. Le vecteur converge vers un vecteur propre de la plus grande valeur propre. Au lieu de cela, l'algorithme QR fonctionne avec une base complète de vecteurs, en utilisant la décomposition QR pour renormaliser (et orthogonaliser). Pour une matrice symétrique A, lors de la convergence, AQ = Q Λ, où Λ est la matrice diagonale des valeurs propres vers lesquelles A a convergé, et où Q est une composition de toutes les transformations de similarité orthogonales nécessaires pour y parvenir. Ainsi, les colonnes de Q sont les vecteurs propres.

Histoire

L'algorithme QR a été précédé par l'algorithme LR, qui utilise la décomposition LU au lieu de la décomposition QR. L'algorithme QR est plus stable, donc l'algorithme LR est rarement utilisé de nos jours. Cependant, cela représente une étape importante dans le développement de l’algorithme QR.

L'algorithme LR a été développé au début des années 1950 par Heinz Rutishauser, qui travaillait à l'époque comme assistant de recherche d'Eduard Stiefel à l'école polytechnique de Zurich. Stiefel a suggéré à Rutishauser d'utiliser la séquence de moments yT
0A^kx₀, k = 0, 1, ... (où x₀ et y₀ sont des vecteurs arbitraires) pour trouver les valeurs propres de A. Rutishauser a pris un algorithme d'Alexander Aitken pour cette tâche et l'a développé en algorithme quotient-différence ou algorithme qd. Après avoir organisé le calcul sous une forme appropriée, il a découvert que l'algorithme qd est en fait l'itération A_k = L_kU_k (décomposition LU), A_{k +1} = U_kL_k, appliquée sur une matrice tridiagonale, d'où découle l'algorithme LR.

Autres variantes

Une variante de l'algorithme QR, l'algorithme de Golub-Kahan-Reinsch, commence par réduire une matrice générale en une matrice bidiagonale. Cette variante de l'algorithme QR pour le calcul de valeur singulières a été d'abord décrit par Golub et Kahan 1965. La routine DBDSQR de la bibliothèque LAPACK implémente cette méthode itérative, avec certaines modifications pour couvrir le cas om les valeurs singulières sont très petites (Demmel et Kahan 1990). Après une première étape utilisant des transformations de Householder et, si approprié, une décomposition QR, on obtient la routine DGESVD pour la décomposition en valeurs singulières. L'algorithme QR peut aussi être implémenté en dimensions infinies avec des résultats de convergence équivalents

Références

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

Sources

• Demmel et Kahan, « Accurate singular values of bidiagonal matrices », SIAM Journal on Scientific and Statistical Computing, vol. 11, n^o 5,‎ 1990, p. 873–912 (DOI 10.1137/0911052, CiteSeer^x 10.1.1.48.3740)
• Golub et Kahan, « Calculating the singular values and pseudo-inverse of a matrix », Journal of the Society for Industrial and Applied Mathematics, Series B: Numerical Analysis, vol. 2, n^o 2,‎ 1965, p. 205–224 (DOI 10.1137/0702016, JSTOR 2949777, Bibcode 1965SJNA....2..205G)

Liens externes

• Eigenvalue problem at PlanetMath.
• Notes on orthogonal bases and the workings of the QR algorithm by Peter J. Olver
• Module for the QR Method