Algèbre

Les décompositions LU et QR

Un ordinateur ne résout jamais un système en calculant un inverse. Il décompose la matrice en un produit de matrices simples : triangulaires (LU: Lower et Upper) ou orthogonale triangulaire (QR). Ces deux décompositions sont le cœur de tout le calcul scientifique. Cinq planches pour les construire pas à pas, et pour voir pourquoi l'ordre des opérations et le choix de la méthode décident de la précision obtenue.

Vue d'ensemble

Résoudre Ax = b est facile quand A est triangulaire : on trouve la dernière inconnue, puis on remonte (substitution). C'est facile aussi quand A est orthogonale (QᵀQ = I) : x = Qᵀb. L'idée est donc d'écrire A comme un produit de telles matrices :

Décomposition	Forme	Fabriquée par	Sert à	Coût pour n × n
LU	PA = LU, L triangulaire inférieure à diagonale 1, U triangulaire supérieure	l'élimination de Gauss, avec échanges de lignes P	systèmes carrés, déterminant, plusieurs seconds membres	≈ 2n³/3
QR	A = QR, Q orthogonale, R triangulaire supérieure	Gram-Schmidt, réflexions de Householder, rotations de Givens	moindres carrés, valeurs propres, bases orthonormées	≈ 4n³/3

L'élimination de Gauss sur une matrice 8 × 8 : chaque étape annule une colonne sous la diagonale. Les multiplicateurs utilisés ne sont pas jetés : ils remplissent L, à gauche. À la fin, A = LU.

LU : l'élimination de Gauss mise en mémoire

Pourquoi garder la trace de l'élimination plutôt que de la refaire ?

Chaque étape de Gauss retranche à la ligne i un multiple m_ik de la ligne k pour annuler un coefficient. En rangeant les m_ik sous la diagonale d'une matrice L (avec des 1 sur la diagonale), on obtient exactement A = LU, où U est la matrice triangulaire finale. Une fois L et U connus, chaque nouveau second membre b se traite en deux substitutions : Ly = b, puis Ux = y.

matrice

À gauche L (multiplicateurs en orange, au fur et à mesure), à droite la matrice en cours d'élimination, qui devient U. La case annulée à l'étape courante est entourée de rouge, la ligne pivot de bleu.

Étape

‖LU − A‖ / ‖A‖ à la fin

det A = produit des pivots

Coût pour k seconds membres (n = taille)

n k seconds membres

Factoriser une fois, résoudre souvent. L'élimination coûte environ 2n³/3 opérations, une substitution seulement 2n². Pour n = 100 et 100 seconds membres, refaire Gauss à chaque fois coûte environ 69 millions d'opérations ; factoriser une fois puis substituer, environ 2,7 millions, soit 25 fois moins. Le compteur ci-dessus a été vérifié en comptant réellement les opérations pour de petites tailles.

Le pivot : pourquoi on échange des lignes

L'élimination sans échange de lignes est exacte en théorie. En machine, elle peut tout perdre.

Prenons A = [[ε, 1], [1, 1]] avec ε minuscule. La solution de Ax = (1 + ε, 2) est exactement x = (1, 1). Sans échange, le multiplicateur vaut 1/ε, gigantesque, et les arrondis de la machine (environ 10⁻¹⁶ en relatif) sont amplifiés d'autant. Avec le pivot partiel, on choisit à chaque étape comme pivot le plus grand coefficient de la colonne : tous les multiplicateurs sont alors au plus égaux à 1 en valeur absolue. On obtient PA = LU, où P est une matrice de permutation.

ε = 10^−k, k

Erreur sur la solution (échelle logarithmique) en fonction de k : sans pivot (rouge), avec pivot partiel (vert). Le point noir marque le k choisi.

Sans pivot : L, U

Sans pivot : x calculé · erreur

Avec pivot : x calculé · erreur

Le pivot partiel n'est pas infaillible

Le danger vient de la croissance des coefficients pendant l'élimination. Avec pivot partiel, ils peuvent au pire doubler à chaque étape. Wilkinson a exhibé la matrice qui atteint ce pire cas : des 1 sur la diagonale et dans la dernière colonne, des −1 sous la diagonale. Sa croissance vaut 2ⁿ⁻¹.

Facteur de croissance (log₂) pour la matrice de Wilkinson (rouge) et maximum observé sur 200 matrices au hasard de même taille (bleu).

Croissance mesurée · 2ⁿ⁻¹

Erreur sur x = (1, …, 1)

Croissance maximale sur 200 matrices au hasard

Un paradoxe célèbre. En théorie, le pivot partiel peut échouer : avec la matrice de Wilkinson, la solution est entièrement fausse dès n = 60, quand la croissance 2⁵⁹ dépasse la précision de la machine. En pratique, sur les matrices rencontrées, la croissance reste modeste (une dizaine pour n = 30, au plus une vingtaine jusqu'à n = 64 sur les matrices tirées ici), et c'est pourquoi tous les logiciels l'utilisent. Comprendre pourquoi les mauvais cas sont si rares reste un problème de recherche.

QR par Gram-Schmidt : orthonormaliser les colonnes

Comment transformer les colonnes de A en une base orthonormée sans perdre l'information ?

On prend les colonnes a₁, a₂, a₃… une par une. On normalise a₁ pour obtenir q₁. On retire à a₂ sa projection sur q₁, puis on normalise : q₂. On retire à a₃ ses projections sur q₁ et q₂, et ainsi de suite. Les coefficients des projections et les normes forment une matrice triangulaire R, et A = QR. C'est le procédé de Gram-Schmidt.

Les vecteurs a₁, a₂, a₃ (gris) et la base q₁, q₂, q₃ qui se construit (couleurs). Les projections retirées apparaissent en pointillés. Faites glisser pour tourner la vue.

Étape

R (triangulaire)

‖QR − A‖ · ‖QᵀQ − I‖

Le même procédé, trois façons de programmer

Gram-Schmidt classique calcule toutes les projections de a_j avec les colonnes d'origine. Gram-Schmidt modifié les retire une à une du vecteur déjà corrigé : c'est mathématiquement identique, mais plus robuste numériquement. Householder (planche 4) procède tout autrement, par réflexions. On les compare sur les matrices de Hilbert hᵢⱼ = 1/(i + j − 1), célèbres pour leur mauvais conditionnement κ.

taille n

Perte d'orthogonalité ‖QᵀQ − I‖ en fonction du conditionnement κ (échelles logarithmiques), pour n = 2 à 12 : classique (rouge), modifié (orange), Householder (vert). Les pointillés gris montrent les pentes κ²·ε et κ·ε.

κ(Hₙ) (valeurs singulières par Jacobi)

Gram-Schmidt classique

Gram-Schmidt modifié

Householder

Ce qu'on mesure. Pour H₁₂ (κ ≈ 1,8·10¹⁶), le Gram-Schmidt classique rend des « q » qui ne sont plus du tout orthogonaux (perte ≈ 1), le modifié perd environ 0,2, et Householder reste à 10⁻¹⁵. La perte du classique croît comme κ², celle du modifié comme κ, celle de Householder pas du tout. Même formule, même résultat exact, mais pas du tout la même précision en machine.

Householder et Givens : annuler par réflexions ou par rotations

Comment faire apparaître des zéros sous la diagonale avec des matrices orthogonales, qui ne déforment rien ?

Une réflexion de Householder H = I − 2vvᵀ (avec ‖v‖ = 1) renvoie un vecteur a sur ‖a‖e₁ par symétrie par rapport à un miroir bien choisi : elle annule d'un coup toute une colonne sous la diagonale. Une rotation de Givens tourne dans le plan de deux coordonnées pour annuler un seul coefficient. Toutes deux conservent les longueurs, donc les erreurs d'arrondi ne sont jamais amplifiées.

transformation

Faites glisser la pointe du vecteur a (bleu). Householder : le miroir (pointillés rouges) envoie a sur (‖a‖, 0). Givens : une rotation d'angle −θ l'amène au même endroit.

Les zéros apparaissent colonne par colonne avec Householder (4 étapes), coefficient par coefficient avec Givens (autant d'étapes que de zéros à créer, ici 14). Rouge : coefficients qui viennent de changer ; gris : zéros.

Étapes

‖QR − A‖ / ‖A‖ · ‖QᵀQ − I‖

Longueur de a avant · après

Lequel choisir ? Householder est le plus rapide pour une matrice pleine : une réflexion par colonne. Givens est préféré quand la matrice est déjà presque triangulaire (creuse, ou de Hessenberg) : on n'annule que les quelques coefficients nécessaires, et les rotations portent sur deux lignes à la fois, ce qui se parallélise bien.

Signature : moindres carrés et valeurs propres

QR résout bien mieux les problèmes de moindres carrés que les équations normales, et, appliquée en boucle, elle trouve les valeurs propres.

On cherche le polynôme de degré d le plus proche, au sens des moindres carrés, de 40 points. La méthode des équations normales résout AᵀAx = Aᵀb : simple, mais le conditionnement de AᵀA est κ(A)², et l'on perd deux fois plus de chiffres. La méthode QR écrit A = QR et résout Rx = Qᵀb, en ne perdant que κ(A). Pour mesurer les erreurs, les données sont tirées exactement d'un polynôme connu, de coefficients 1, 1/2, 1/3, …

degré d

Erreur sur les coefficients (échelle logarithmique) en fonction du degré : équations normales (rouge), QR (vert). En gris, κ(A)·ε et κ(A)²·ε.

κ(A) (Jacobi)

Erreur des équations normales

Erreur de QR (Householder)

Rapport

On factorise A₀ = Q₀R₀, puis on multiplie dans l'autre ordre : A₁ = R₀Q₀ = Q₀ᵀA₀Q₀. Cette nouvelle matrice est semblable à A₀, donc elle a les mêmes valeurs propres. En recommençant, les coefficients sous la diagonale s'éteignent et la diagonale converge vers les valeurs propres, rangées par ordre décroissant de module. La matrice symétrique ci-dessous cache les valeurs propres 5, 4, 3, 2 et 1.

|aᵢⱼ| en couleur (échelle logarithmique, du clair 10⁻¹⁶ au foncé), diagonale encadrée. Les coefficients hors diagonale disparaissent au fil des itérations.

log₁₀ |a_i+1,i| en fonction de l'itération, pour i = 1 à 4. Les pointillés ont la pente prévue log₁₀(λ_i+1/λ_i) : 4/5, 3/4, 2/3, 1/2.

Itération

Diagonale

Taux mesurés |a_i+1,i(k+1)/a_i+1,i(k)| · prévus λ_i+1/λ_i

L'algorithme du siècle. Proposé en 1961 par John Francis et Vera Kublanovskaïa, l'algorithme QR est le moyen standard de calculer toutes les valeurs propres d'une matrice. Sous sa forme simple montrée ici, il converge au rythme |λ_i+1/λ_i|, lent si deux valeurs propres sont proches. Les versions réelles ajoutent des « décalages » qui rendent la convergence quadratique, voire cubique. Côté moindres carrés, le rapport d'erreur entre équations normales et QR atteint ici plusieurs ordres de grandeur dès le degré 10.

Une brève histoire

1809–1810Gauss publie la méthode des moindres carrés (Theoria motus), puis résout les équations normales par élimination pour calculer l'orbite de l'astéroïde Pallas.

1907Erhard Schmidt publie le procédé d'orthonormalisation, déjà utilisé par Gram (1883) et Laplace.

1947Von Neumann et Goldstine analysent les erreurs d'arrondi de l'élimination de Gauss sur les premiers ordinateurs.

1948Turing introduit la décomposition LU et la notion de conditionnement.

1958Householder propose ses réflexions pour la factorisation QR. Givens avait publié ses rotations en 1954.

1961Wilkinson analyse la croissance des coefficients avec pivot partiel. La même année, Francis et Kublanovskaïa inventent l'algorithme QR.

En résumé

LU

L'élimination de Gauss, avec ses multiplicateurs gardés dans L. On factorise une fois (2n³/3) et l'on résout ensuite chaque système en 2n². Le pivot partiel est indispensable en machine.

QR

Orthogonale fois triangulaire. Gram-Schmidt est la version des manuels, Householder et Givens celles des logiciels, car les matrices orthogonales n'amplifient pas les erreurs.

Usages

Les moindres carrés par QR perdent κ chiffres au lieu de κ². Itérer QR, puis RQ, fait converger la matrice vers ses valeurs propres.

Pour aller plus loin

Cholesky. Pour une matrice symétrique définie positive, A = CCᵀ avec C triangulaire : deux fois moins coûteux que LU, et stable sans pivot.

Pivot total, pivot de colonne. Le pivot total (lignes et colonnes) borne mieux la croissance, mais coûte plus cher. La factorisation QR avec pivot de colonnes révèle le rang numérique d'une matrice.

Matrices creuses. Pour les très grandes matrices issues des équations aux dérivées partielles, l'ordre d'élimination décide du « remplissage » de L et U ; on combine alors décompositions et méthodes itératives.

Liens. Systèmes linéaires (VIII.14), matrices et valeurs singulières, espaces euclidiens (VIII.22) pour l'orthogonalité, et la page sur les valeurs singulières pour le conditionnement.