Décompositions matricielles et calcul numérique

Objectifs du chapitre

À la fin de ce chapitre, vous serez capable de:

• calculer la décomposition $A=LU$ d'une matrice carrée, l'utiliser pour résoudre plusieurs systèmes de même matrice et pour calculer un déterminant, et expliquer pourquoi le pivotage partiel est indispensable en arithmétique flottante;
• calculer la décomposition de Cholesky $A=L{L}^{\mathsf{T}}$ d'une matrice symétrique définie positive;
• obtenir la décomposition $A=QR$ par le procédé de Gram–Schmidt ou par des réflexions de Householder, et l'appliquer aux moindres carrés et au calcul des valeurs propres;
• énoncer et démontrer le théorème de décomposition en valeurs singulières $A=U\Sigma V^{\mathsf{T}}$, l'interpréter géométriquement et en lire le rang, les quatre sous-espaces fondamentaux, les normes et la pseudo-inverse;
• définir le conditionnement $\kappa (A)$, démontrer l'inégalité $\Vert \delta \mathbf{x}\Vert /\Vert \mathbf{x}\Vert \le \kappa (A)\Vert \delta \mathbf{b}\Vert /\Vert \mathbf{b}\Vert$ et en tirer les bonnes pratiques du calcul numérique;
• décrire les applications de la SVD tronquée (compression, analyse en composantes principales) et les méthodes itératives de Jacobi et de Gauss–Seidel.

Pourquoi factoriser une matrice?

Les chapitres précédents ont construit les objets de l'algèbre linéaire: systèmes, matrices, déterminants, espaces vectoriels, orthogonalité, valeurs propres. Ce dernier chapitre change de point de vue et pose la question de l'ingénieur: comment calcule-t-on tout cela, vite et sûrement, sur un ordinateur? Quelques situations typiques.

• Des milliers de systèmes de même matrice. Dans une simulation par éléments finis, la matrice de rigidité $K$ d'une structure ne change pas d'un pas de temps à l'autre; seul le second membre (les charges) change. Résoudre $K\mathbf{x}={\mathbf{b}}_k$ pour $k=1,\dots ,10000$ en refaisant chaque fois l'élimination de Gauss serait un gaspillage: on veut factoriser une fois, résoudre souvent.
• Ajuster des données bruitées. Un géomètre mesure des distances et des angles entre des dizaines de points; les inconnues (les coordonnées) sont bien moins nombreuses que les mesures. Le système est surdéterminé et il faut le résoudre au sens des moindres carrés (chapitre 7), sans que les erreurs de mesure soient amplifiées par le calcul.
• Compresser une image. Une photographie de $1000\times 1000$ pixels est une matrice d'un million de coefficients. Peut-on la remplacer par une matrice de rang $50$, stockée avec dix fois moins de nombres, sans perte visible?
• Trouver les directions principales d'un nuage de points. Des capteurs enregistrent cent grandeurs à chaque instant. Ces cent colonnes sont-elles vraiment indépendantes, ou le phénomène se décrit-il avec trois combinaisons bien choisies?

À chacune de ces questions répond une factorisation: écrire $A$ comme un produit de matrices plus simples (triangulaires, orthogonales, diagonales), pour lesquelles toutes les opérations sont faciles. Les décompositions $LU$, $QR$ et en valeurs singulières (SVD) sont les outils de base de l'algèbre linéaire numérique; elles constituent le cœur de la bibliothèque LAPACK, que NumPy, MATLAB, R ou Julia appellent en coulisses. Comprendre ce qu'elles font, ce qu'elles coûtent et quand elles échouent est indispensable pour utiliser ces logiciels intelligemment.

La décomposition LU

De l'élimination de Gauss à A = LU

Revenons à l'élimination de Gauss du chapitre 1 sur une matrice carrée $A\in {\mathcal{M}}_n(\mathbb{R})$, en supposant que l'on n'ait jamais besoin d'échanger deux lignes. Pour annuler le coefficient $a_{ij}$ sous le pivot $a_{jj}$, on effectue $L_i\leftarrow L_i-{\ell }_{ij}{L}_j$ avec le multiplicateur ${\ell }_{ij}=a_{ij}/a_{jj}$ (calculé sur la matrice courante). Au chapitre 2, nous avons vu que cette opération revient à multiplier à gauche par la matrice élémentaire $E_{ij}=I-{\ell }_{ij}{\mathbf{e}}_i{\mathbf{e}}_j^{\mathsf{T}}$, triangulaire inférieure à diagonale de $1$. À la fin de l'élimination, on a obtenu une matrice triangulaire supérieure $U$:

Démonstration. Chaque $E_{ij}$ est inversible d'inverse $E_{ij}^{-1}=I+{\ell }_{ij}{\mathbf{e}}_i{\mathbf{e}}_j^{\mathsf{T}}$ (on vérifie que $(I-\ell {\mathbf{e}}_i{\mathbf{e}}_j^{\mathsf{T}})(I+\ell {\mathbf{e}}_i{\mathbf{e}}_j^{\mathsf{T}})=I-{\ell }^2{\mathbf{e}}_i({\mathbf{e}}_j^{\mathsf{T}}{\mathbf{e}}_i){\mathbf{e}}_j^{\mathsf{T}}=I$ car ${\mathbf{e}}_j^{\mathsf{T}}{\mathbf{e}}_i=0$ pour $i\ne j$). Donc

le produit étant pris dans l'ordre de l'élimination (colonne $1$ d'abord). Développons ce produit: un terme croisé $({\ell }_{ij}{\mathbf{e}}_i{\mathbf{e}}_j^{\mathsf{T}})({\ell }_{kl}{\mathbf{e}}_k{\mathbf{e}}_l^{\mathsf{T}})$ contient le scalaire ${\mathbf{e}}_j^{\mathsf{T}}{\mathbf{e}}_k$, nul dès que $j\ne k$; or dans l'ordre choisi, le facteur $(i,j)$ précède le facteur $(k,l)$ seulement si $j\le l<k$, de sorte que $j\ne k$ toujours. Tous les termes croisés s'annulent et

qui est bien la matrice triangulaire inférieure unitaire contenant les multiplicateurs. Pour l'unicité, si $A=L_1{U}_1=L_2{U}_2$ avec $A$ inversible, alors $U_1,U_2$ sont inversibles (leurs déterminants sont non nuls) et $L_2^{-1}{L}_1=U_2{U}_1^{-1}$. Le membre de gauche est triangulaire inférieur unitaire (produit et inverse de telles matrices le sont, chapitre 2), le membre de droite triangulaire supérieur: la seule matrice qui soit les deux à la fois est $I$. Donc $L_1=L_2$ et $U_1=U_2$. $\square$

Résoudre A x = b: substitution avant, substitution arrière

Une fois $A=LU$ connue, le système $A\mathbf{x}=\mathbf{b}$ s'écrit $L(U\mathbf{x})=\mathbf{b}$ et se résout en deux systèmes triangulaires:

Le premier se résout de haut en bas ($y_1=b_1$, puis $y_2=b_2-{\ell }_{21}{y}_1$, etc.), le second de bas en haut.

Coût. Comptons les opérations flottantes (additions et multiplications, flops). L'élimination de la colonne $j$ modifie un bloc de $(n-j)$ lignes et $(n-j)$ colonnes, à raison de deux flops par coefficient: au total ${\sum }_{j=1}^{n-1}2(n-j{)}^2\approx \frac{2n^3}{3}$ flops pour la factorisation. Chaque substitution triangulaire coûte environ $n^2$ flops, soit $2n^2$ par second membre. Pour $n=1000$, la factorisation demande environ $6,7\cdot 10^8$ flops et chaque résolution $2\cdot 10^6$: on peut résoudre plus de trois cents systèmes pour le prix d'une factorisation. Calculer $A^{-1}$ puis multiplier $A^{-1}\mathbf{b}$ coûterait trois fois plus cher ($2n^3$ flops pour l'inverse) et, nous le verrons, serait moins précis.

Le pivotage partiel

L'hypothèse «aucun échange de lignes» n'est pas seulement une commodité théorique: la matrice $\left(\begin{array}{cc}0& 1\\ 1& 1\end{array}\right)$, inversible, n'a pas de décomposition LU (le premier pivot est nul). Plus grave, en arithmétique flottante, un pivot petit est presque aussi dangereux qu'un pivot nul.

La parade est le pivotage partiel: à chaque étape, on choisit comme pivot le coefficient de plus grande valeur absolue dans la colonne courante (sur et sous la diagonale) et l'on échange les lignes en conséquence. Tous les multiplicateurs vérifient alors $|{\ell }_{ij}|\le 1$, ce qui empêche l'amplification catastrophique ci-dessus. Les échanges de lignes se résument en une matrice de permutation $P$ (chapitre 2).

Nous admettons ce théorème; l'idée est simple. À chaque étape, la colonne courante (sous la diagonale incluse) n'est pas entièrement nulle, sinon la matrice serait singulière: le pivot de plus grande valeur absolue est donc non nul. On montre ensuite qu'un échange de lignes effectué à l'étape $k$ peut être «remonté» avant toutes les opérations d'élimination précédentes, à condition d'échanger aussi les lignes correspondantes des multiplicateurs déjà stockés dans $L$; toutes les permutations se regroupent ainsi dans $P$, et l'on est ramené au théorème 12.1 pour la matrice $PA$. Le système $A\mathbf{x}=\mathbf{b}$ devient $LU\mathbf{x}=P\mathbf{b}$: on permute le second membre, puis on substitue comme en (12.1). C'est exactement ce que fait scipy.linalg.lu_factor ou l'opérateur \ de MATLAB pour une matrice carrée quelconque.

Déterminant via LU

Le déterminant d'un produit est le produit des déterminants (chapitre 3), celui d'une matrice triangulaire est le produit de ses coefficients diagonaux, et $\det P=\pm 1$ selon la parité du nombre d'échanges. De $PA=LU$ on tire donc

Pour la matrice de l'exemple 12.1, $\det A=2\cdot 1\cdot 2=4$ (aucun échange). C'est ainsi que numpy.linalg.det calcule un déterminant: en $2n^3/3$ flops, et non par le développement selon une ligne, dont le coût en $n!$ est prohibitif dès $n=20$.

La décomposition de Cholesky

Lorsque $A$ est symétrique définie positive — cas des matrices de rigidité en mécanique des structures, des matrices $A^{\mathsf{T}}A$ des moindres carrés, des matrices de covariance en statistique — on peut faire deux fois mieux.

Démonstration. Au chapitre 11, nous avons vu qu'une matrice symétrique est définie positive si et seulement si tous ses mineurs principaux dominants ${\Delta }_k=\det A_k$ ($A_k$ le bloc $k\times k$ en haut à gauche) sont strictement positifs. L'élimination de Gauss sans échange sur $A$ transforme $A_k$ en le bloc $U_k$, sans changer son déterminant: ${\Delta }_k=u_{11}\cdots u_{kk}$, d'où le $k$-ième pivot $u_{kk}={\Delta }_k/{\Delta }_{k-1}>0$. Aucun pivot n'est nul et le théorème 12.1 donne $A=LU$ de manière unique. Notons $D=\mathrm{diag}(u_{11},\dots ,u_{nn})$ et $\tilde{U}=D^{-1}U$, triangulaire supérieure unitaire: $A=LD\tilde{U}$. La symétrie de $A$ donne $A=A^{\mathsf{T}}={\tilde{U}}^{\mathsf{T}}D{L}^{\mathsf{T}}={\tilde{U}}^{\mathsf{T}}(D{L}^{\mathsf{T}})$, qui est une autre décomposition LU de $A$ (facteur inférieur unitaire ${\tilde{U}}^{\mathsf{T}}$, facteur supérieur $D{L}^{\mathsf{T}}$). Par unicité, ${\tilde{U}}^{\mathsf{T}}=L$, donc $A=LD{L}^{\mathsf{T}}$. Comme les $u_{kk}$ sont positifs, $D^{1/2}=\mathrm{diag}(\sqrt{u_{kk}})$ existe et $A=(L{D}^{1/2})(L{D}^{1/2}{)}^{\mathsf{T}}$: le facteur $L{D}^{1/2}$ est triangulaire inférieur à diagonale $\sqrt{u_{kk}}>0$. L'unicité découle de celle de la décomposition LU: si $A=L_1{L}_1^{\mathsf{T}}$, alors $A=(L_1{D}_1^{-1})(D_1{L}_1^{\mathsf{T}})$ avec $D_1$ la diagonale de $L_1$, ce qui force $L_1{D}_1^{-1}=L$ et $D_1^2=D$. $\square$

Algorithme. On identifie les coefficients de $A=L{L}^{\mathsf{T}}$ colonne par colonne, ce qui donne, pour $j=1,\dots ,n$,

Le coût est d'environ $n^3/3$ flops, la moitié de la décomposition LU, et aucun pivotage n'est nécessaire (la stabilité est automatique: on montre que ${\ell }_{ik}^2\le a_{ii}$). Si une racine carrée d'un nombre négatif apparaît, c'est que $A$ n'était pas définie positive: l'algorithme de Cholesky est d'ailleurs le test numérique standard de positivité.

La décomposition QR

Gram–Schmidt vu comme une factorisation

Le procédé de Gram–Schmidt (chapitre 7) transforme une famille libre $({\mathbf{a}}_1,\dots ,{\mathbf{a}}_n)$ de ${\mathbb{R}}^m$ en une famille orthonormée $({\mathbf{q}}_1,\dots ,{\mathbf{q}}_n)$ qui engendre les mêmes sous-espaces emboîtés: $\mathrm{Vect}({\mathbf{q}}_1,\dots ,{\mathbf{q}}_j)=\mathrm{Vect}({\mathbf{a}}_1,\dots ,{\mathbf{a}}_j)$ pour tout $j$. Rangeons les ${\mathbf{a}}_j$ en colonnes d'une matrice $A\in {\mathcal{M}}_{m\times n}(\mathbb{R})$ et les ${\mathbf{q}}_j$ en colonnes de $Q$.

Démonstration. Le procédé de Gram–Schmidt construit, pour $j=1,\dots ,n$,

avec ${\mathbf{w}}_j\ne 0$ parce que ${\mathbf{a}}_j\notin \mathrm{Vect}({\mathbf{a}}_1,\dots ,{\mathbf{a}}_{j-1})=\mathrm{Vect}({\mathbf{q}}_1,\dots ,{\mathbf{q}}_{j-1})$. En posant $r_{ij}=⟨{\mathbf{q}}_i,{\mathbf{a}}_j⟩$ pour $i<j$, $r_{jj}=\Vert {\mathbf{w}}_j\Vert >0$ et $r_{ij}=0$ pour $i>j$, la relation ci-dessus se réécrit

ce qui est exactement la $j$-ième colonne de l'égalité matricielle $A=QR$. $\square$

Les réflexions de Householder

Le procédé de Gram–Schmidt est parfait sur le papier, mais en arithmétique flottante il perd progressivement l'orthogonalité: si deux colonnes sont presque parallèles, la soustraction ${\mathbf{a}}_j-\sum ⟨{\mathbf{q}}_i,{\mathbf{a}}_j⟩{\mathbf{q}}_i$ est une différence de nombres proches, et les erreurs d'arrondi dominent le petit reste ${\mathbf{w}}_j$. Les bibliothèques numériques calculent $QR$ autrement: en appliquant à $A$ une suite de réflexions orthogonales qui annulent, colonne après colonne, les coefficients sous la diagonale, exactement comme l'élimination de Gauss, mais avec des matrices orthogonales à la place des matrices élémentaires.

Démonstration. 1. Posons $P=\mathbf{v}{\mathbf{v}}^{\mathsf{T}}/({\mathbf{v}}^{\mathsf{T}}\mathbf{v})$, la matrice de la projection orthogonale sur la droite $\mathbb{R}\mathbf{v}$ (chapitre 7). Elle est symétrique et $P^2=\mathbf{v}({\mathbf{v}}^{\mathsf{T}}\mathbf{v}){\mathbf{v}}^{\mathsf{T}}/({\mathbf{v}}^{\mathsf{T}}\mathbf{v}{)}^2=P$. Donc $H_{\mathbf{v}}=I-2P$ est symétrique et

ce qui prouve à la fois l'orthogonalité et $H_{\mathbf{v}}^2=I$.

2. Notons $\alpha =\Vert \mathbf{x}\Vert$ et $\mathbf{v}=\mathbf{x}-\alpha {\mathbf{e}}_1\ne 0$. Alors ${\mathbf{v}}^{\mathsf{T}}\mathbf{x}={\alpha }^2-\alpha x_1$ et ${\mathbf{v}}^{\mathsf{T}}\mathbf{v}={\alpha }^2-2\alpha x_1+{\alpha }^2=2({\alpha }^2-\alpha x_1)=2{\mathbf{v}}^{\mathsf{T}}\mathbf{x}$. Par conséquent

Géométriquement, $\mathbf{x}$ et $\alpha {\mathbf{e}}_1$ ont la même longueur; la réflexion par rapport à l'hyperplan médiateur de ces deux vecteurs (orthogonal à leur différence $\mathbf{v}$) envoie l'un sur l'autre. Par exemple, pour $\mathbf{x}=(3,4)$, $\alpha =5$, $\mathbf{v}=(-2,4)$, ${\mathbf{v}}^{\mathsf{T}}\mathbf{v}=20$ et

L'algorithme de Householder. On choisit $H_1$ qui envoie la première colonne de $A$ sur $\Vert {\mathbf{a}}_1\Vert {\mathbf{e}}_1$: la première colonne de $H_{1}A$ est nulle sous la diagonale. On recommence avec une réflexion $H_2$ agissant sur les lignes $2$ à $m$ (bloc-diagonale avec $1$ en position $(1,1)$) pour annuler la deuxième colonne sous la diagonale, sans toucher à la première, et ainsi de suite. Après $n$ étapes, $H_n\cdots H_2{H}_{1}A=\left(\begin{array}{c}R\\ 0\end{array}\right)$, d'où $A=Q_{\text{c}}\left(\begin{array}{c}R\\ 0\end{array}\right)$ avec $Q_{\text{c}}=H_1{H}_2\cdots H_n$ orthogonale $m\times m$ (chaque $H_k$ est son propre inverse). Les $n$ premières colonnes de $Q_{\text{c}}$ forment le $Q$ du théorème 12.4 (au signe des colonnes près).

Application aux moindres carrés

Au chapitre 7, la solution au sens des moindres carrés de $A\mathbf{x}=\mathbf{b}$ (système surdéterminé, $m>n$) était obtenue par les équations normales $A^{\mathsf{T}}A\hat{\mathbf{x}}=A^{\mathsf{T}}\mathbf{b}$. La décomposition QR en donne une version plus stable.

Démonstration. Comme $A$ est de rang $n$, $A^{\mathsf{T}}A$ est inversible et les équations normales ont une solution unique (chapitre 7). Remplaçons $A$ par $QR$: $A^{\mathsf{T}}A=R^{\mathsf{T}}{Q}^{\mathsf{T}}QR=R^{\mathsf{T}}R$ puisque $Q^{\mathsf{T}}Q=I_n$, et $A^{\mathsf{T}}\mathbf{b}=R^{\mathsf{T}}{Q}^{\mathsf{T}}\mathbf{b}$. Les équations normales s'écrivent donc $R^{\mathsf{T}}R\hat{\mathbf{x}}=R^{\mathsf{T}}{Q}^{\mathsf{T}}\mathbf{b}$; comme $R$ est inversible (diagonale strictement positive), $R^{\mathsf{T}}$ l'est aussi et l'on peut simplifier: $R\hat{\mathbf{x}}=Q^{\mathsf{T}}\mathbf{b}$. $\square$

L'algorithme QR pour les valeurs propres

La décomposition QR a une seconde application, surprenante: elle est au cœur de l'algorithme standard de calcul des valeurs propres, celui de numpy.linalg.eig. Le polynôme caractéristique, si commode au chapitre 9, est inutilisable numériquement: ses coefficients sont extrêmement sensibles aux arrondis et ses racines le sont encore plus.

Principe. Partant de $A_0=A$, on factorise $A_k=Q_k{R}_k$ puis on multiplie les facteurs dans l'ordre inverse:

Comme $R_k=Q_k^{\mathsf{T}}{A}_k$, on a $A_{k+1}=Q_k^{\mathsf{T}}{A}_k{Q}_k$: toutes les matrices $A_k$ sont orthogonalement semblables à $A$ et ont les mêmes valeurs propres (chapitre 6). Le fait remarquable, que nous admettons, est le suivant: si $A$ est symétrique à valeurs propres de valeurs absolues distinctes, $A_k$ converge vers la matrice diagonale des valeurs propres, rangées par valeur absolue décroissante, et les produits $Q_0{Q}_1\cdots Q_k$ convergent vers la matrice orthogonale des vecteurs propres. Pour une matrice quelconque, $A_k$ tend vers une forme triangulaire supérieure (forme de Schur), dont la diagonale porte les valeurs propres.

La décomposition en valeurs singulières

Les décompositions LU et QR sont des outils de calcul. La décomposition en valeurs singulières (SVD, singular value decomposition) est davantage: elle révèle la structure géométrique complète d'une matrice quelconque, même rectangulaire, même de rang déficient, là où la diagonalisation exige une matrice carrée et diagonalisable.

Le théorème

Démonstration. La matrice $A^{\mathsf{T}}A\in {\mathcal{M}}_n(\mathbb{R})$ est symétrique, et semi-définie positive: pour tout $\mathbf{x}$, ${\mathbf{x}}^{\mathsf{T}}{A}^{\mathsf{T}}A\mathbf{x}=\Vert A\mathbf{x}{\Vert }^2\ge 0$. Par le théorème spectral (chapitre 11), elle possède une base orthonormée de vecteurs propres ${\mathbf{v}}_1,\dots ,{\mathbf{v}}_n$ associés à des valeurs propres ${\lambda }_1\ge \cdots \ge {\lambda }_n\ge 0$ (positives par la remarque précédente appliquée à $\mathbf{x}={\mathbf{v}}_i$: ${\lambda }_i=\Vert A{\mathbf{v}}_i{\Vert }^2$). Soit $r$ le nombre de valeurs propres strictement positives; posons ${\sigma }_i=\sqrt{{\lambda }_i}$ et, pour $i\le r$,

Ces vecteurs sont orthonormés: $⟨{\mathbf{u}}_i,{\mathbf{u}}_j⟩=\frac{{\mathbf{v}}_i^{\mathsf{T}}{A}^{\mathsf{T}}A{\mathbf{v}}_j}{{\sigma }_i{\sigma }_j}=\frac{{\lambda }_j{\mathbf{v}}_i^{\mathsf{T}}{\mathbf{v}}_j}{{\sigma }_i{\sigma }_j}={\delta }_{ij}$. On complète $({\mathbf{u}}_1,\dots ,{\mathbf{u}}_r)$ en une base orthonormée $({\mathbf{u}}_1,\dots ,{\mathbf{u}}_m)$ de ${\mathbb{R}}^m$ (théorème de la base incomplète puis Gram–Schmidt), rangée dans $U$. Pour $i>r$, $\Vert A{\mathbf{v}}_i{\Vert }^2={\lambda }_i=0$ donc $A{\mathbf{v}}_i=0$. Au total, $A{\mathbf{v}}_i={\sigma }_i{\mathbf{u}}_i$ pour $i\le r$ et $A{\mathbf{v}}_i=0$ sinon, ce qui s'écrit matriciellement $AV=U\Sigma$, et $V$ étant orthogonale, $A=U\Sigma V^{\mathsf{T}}$. $\square$

La démonstration est constructive, mais ce n'est pas ainsi que les logiciels calculent la SVD: former $A^{\mathsf{T}}A$ élève le conditionnement au carré (théorème 12.13 plus bas). L'algorithme de Golub et Kahan (1965) réduit $A$ à une forme bidiagonale par des réflexions de Householder, puis applique une variante de l'algorithme QR.

Interprétation géométrique

Lisons $A=U\Sigma V^{\mathsf{T}}$ de droite à gauche sur un vecteur $\mathbf{x}$: $V^{\mathsf{T}}$ est une isométrie (rotation ou réflexion) qui envoie la base orthonormée $({\mathbf{v}}_i)$ sur la base canonique; $\Sigma$ étire l'axe $i$ d'un facteur ${\sigma }_i$ (et écrase les axes $i>r$); $U$ est une nouvelle isométrie qui envoie la base canonique sur $({\mathbf{u}}_i)$. Toute application linéaire est une rotation, suivie d'un étirement le long d'axes orthogonaux, suivie d'une rotation. En particulier, l'image de la sphère unité de ${\mathbb{R}}^n$ est un ellipsoïde de ${\mathbb{R}}^m$ (éventuellement aplati) dont les demi-axes sont les vecteurs ${\sigma }_i{\mathbf{u}}_i$.

Rang et sous-espaces fondamentaux

Démonstration. Par (12.11), pour tout $\mathbf{x}$, $A\mathbf{x}={\sum }_{i=1}^r{\sigma }_i⟨{\mathbf{v}}_i,\mathbf{x}⟩{\mathbf{u}}_i\in \mathrm{Vect}({\mathbf{u}}_1,\dots ,{\mathbf{u}}_r)$, et réciproquement chaque ${\mathbf{u}}_i=A({\mathbf{v}}_i/{\sigma }_i)$ est dans l'image: $\mathrm{Im}A=\mathrm{Vect}({\mathbf{u}}_1,\dots ,{\mathbf{u}}_r)$, de dimension $r$ puisque les ${\mathbf{u}}_i$ sont orthonormés donc libres. D'autre part, $A{\mathbf{v}}_i=0$ pour $i>r$, donc $\mathrm{Vect}({\mathbf{v}}_{r+1},\dots ,{\mathbf{v}}_n)\subseteq \ker A$; par le théorème du rang (chapitre 5), $\dim \ker A=n-r$, et l'inclusion est une égalité. Les deux autres formules s'obtiennent en appliquant ce qui précède à $A^{\mathsf{T}}=V{\Sigma }^{\mathsf{T}}{U}^{\mathsf{T}}$, qui est une SVD de $A^{\mathsf{T}}$ avec les rôles de $U$ et $V$ échangés. $\square$

Ce théorème donne des bases orthonormées des quatre sous-espaces fondamentaux d'un coup, et rend visible le théorème du chapitre 7 selon lequel $\ker A=(\mathrm{Im}{A}^{\mathsf{T}}{)}^{\perp }$ et $\ker A^{\mathsf{T}}=(\mathrm{Im}A{)}^{\perp }$. Par exemple, pour $A=\left(\begin{array}{cc}1& 2\\ 2& 4\end{array}\right)$, on trouve $A^{\mathsf{T}}A=\left(\begin{array}{cc}5& 10\\ 10& 20\end{array}\right)$, de valeurs propres $25$ et $0$: une seule valeur singulière ${\sigma }_1=5$, donc $\mathrm{rg}A=1$; ${\mathbf{v}}_1=(1,2)/\sqrt{5}$ engendre $\mathrm{Im}{A}^{\mathsf{T}}$ et ${\mathbf{v}}_2=(2,-1)/\sqrt{5}$ engendre $\ker A$; ${\mathbf{u}}_1=A{\mathbf{v}}_1/5=(1,2)/\sqrt{5}$ engendre $\mathrm{Im}A$ et ${\mathbf{u}}_2=(2,-1)/\sqrt{5}$ engendre $\ker A^{\mathsf{T}}$. Numériquement, le rang d'une matrice se calcule en comptant les valeurs singulières supérieures à un seuil (par défaut ${\sigma }_1\cdot \max (m,n)\cdot {\epsilon }_{\text{machine}}$ dans NumPy): l'élimination de Gauss, elle, ne sait pas distinguer un pivot nul d'un pivot de $10^{-15}$.

Normes matricielles

Démonstration. Soit $\mathbf{x}$ unitaire et $\mathbf{y}=V^{\mathsf{T}}\mathbf{x}$, unitaire lui aussi car $V$ est orthogonale. Alors $\Vert A\mathbf{x}{\Vert }^2=\Vert U\Sigma \mathbf{y}{\Vert }^2=\Vert \Sigma \mathbf{y}{\Vert }^2={\sum }_{i=1}^r{\sigma }_i^2{y}_i^2\le {\sigma }_1^2{\sum }_i{y}_i^2={\sigma }_1^2$, avec égalité pour $\mathbf{y}={\mathbf{e}}_1$, c'est-à-dire $\mathbf{x}={\mathbf{v}}_1$. Pour Frobenius, $\Vert A{\Vert }_F^2=\mathrm{tr}(A^{\mathsf{T}}A)=\mathrm{tr}(V{\Sigma }^{\mathsf{T}}\Sigma V^{\mathsf{T}})=\mathrm{tr}({\Sigma }^{\mathsf{T}}\Sigma )=\sum {\sigma }_i^2$, en utilisant $\mathrm{tr}(MN)=\mathrm{tr}(NM)$ (chapitre 9). $\square$

La pseudo-inverse

Démonstration. Posons $\mathbf{y}=V^{\mathsf{T}}\mathbf{x}$ et $\mathbf{c}=U^{\mathsf{T}}\mathbf{b}$. Comme $U$ est orthogonale, $\Vert A\mathbf{x}-\mathbf{b}{\Vert }^2=\Vert U(\Sigma \mathbf{y}-\mathbf{c}){\Vert }^2={\sum }_{i=1}^r({\sigma }_i{y}_i-c_i{)}^2+{\sum }_{i=r+1}^m{c}_i^2$. Le second terme ne dépend pas de $\mathbf{x}$; le premier est nul, donc minimal, si et seulement si $y_i=c_i/{\sigma }_i$ pour $i\le r$, les composantes $y_{r+1},\dots ,y_n$ restant libres. Parmi ces minimiseurs, $\Vert \mathbf{x}\Vert =\Vert \mathbf{y}\Vert$ est minimale lorsque $y_i=0$ pour $i>r$; le vecteur obtenu est exactement $\mathbf{y}={\Sigma }^{+}\mathbf{c}$, soit $\mathbf{x}=V{\Sigma }^{+}{U}^{\mathsf{T}}\mathbf{b}=A^{+}\mathbf{b}$. $\square$

Reprenons $A=\left(\begin{array}{cc}1& 2\\ 2& 4\end{array}\right)$ et $\mathbf{b}=(1,0)$. Le système $A\mathbf{x}=\mathbf{b}$ n'a pas de solution ($\mathbf{b}\notin \mathrm{Im}A$) et, $A$ étant de rang $1$, les solutions au sens des moindres carrés forment une droite entière. La pseudo-inverse est $A^{+}=\frac{1}{{\sigma }_1}{\mathbf{v}}_1{\mathbf{u}}_1^{\mathsf{T}}=\frac{1}{5}\cdot \frac{1}{5}\left(\begin{array}{c}1\\ 2\end{array}\right)\left(\begin{array}{cc}1& 2\end{array}\right)=\frac{1}{25}\left(\begin{array}{cc}1& 2\\ 2& 4\end{array}\right)$, d'où ${\mathbf{x}}^{+}=A^{+}\mathbf{b}=\frac{1}{25}(1,2)$. On vérifie que $A{\mathbf{x}}^{+}=\frac{1}{5}(1,2)$ est bien la projection orthogonale de $\mathbf{b}$ sur $\mathrm{Im}A=\mathbb{R}{\mathbf{u}}_1$, et que ${\mathbf{x}}^{+}$ est orthogonal à $\ker A=\mathbb{R}(2,-1)$, ce qui en fait le minimiseur le plus court. C'est ce que calcule numpy.linalg.lstsq pour un système de rang déficient.

Conditionnement et erreurs d'arrondi

Nombres flottants

Un ordinateur représente les réels en virgule flottante: en double précision (norme IEEE 754, type float de Python), un nombre s'écrit $\pm 1,m\times 2^e$ avec une mantisse $m$ de 52 bits. Entre deux flottants consécutifs, l'écart relatif est l'epsilon machine $\epsilon \approx 2,2\cdot 10^{-16}$: chaque opération arithmétique élémentaire est effectuée avec une erreur relative au plus $\epsilon /2$, soit environ 16 chiffres significatifs justes. Cela paraît beaucoup, mais nous avons vu avec le pivot $10^{-17}$ qu'un algorithme mal conçu peut perdre tous ces chiffres en une seule soustraction. La question centrale est la suivante: en quoi la précision de la solution dépend-elle de la précision des données?

Le conditionnement d'une matrice

Démonstration. Si $A=U\Sigma V^{\mathsf{T}}$ est inversible, alors $A^{-1}=V{\Sigma }^{-1}{U}^{\mathsf{T}}$, qui est une SVD de $A^{-1}$ (matrices orthogonales $V$ et $U$, matrice diagonale ${\Sigma }^{-1}$ de coefficients $1/{\sigma }_n\ge \cdots \ge 1/{\sigma }_1$, à réordonner). Par la proposition 12.9, $\Vert A{\Vert }_2={\sigma }_1$ et $\Vert A^{-1}{\Vert }_2=1/{\sigma }_n$. $\square$

Géométriquement, ${\kappa }_2(A)$ est le rapport du grand au petit demi-axe de l'ellipsoïde image de la sphère unité (explorateur ci-dessus): une matrice mal conditionnée aplatit la sphère en une galette, et l'on ne peut plus distinguer, dans l'image, des directions qui étaient très différentes au départ. Le théorème suivant transforme cette image en une borne d'erreur.

Démonstration. Par soustraction des deux systèmes, $A\delta \mathbf{x}=\delta \mathbf{b}$, donc $\delta \mathbf{x}=A^{-1}\delta \mathbf{b}$ et $\Vert \delta \mathbf{x}\Vert \le \Vert A^{-1}\Vert \Vert \delta \mathbf{b}\Vert$. D'autre part $\Vert \mathbf{b}\Vert =\Vert A\mathbf{x}\Vert \le \Vert A\Vert \Vert \mathbf{x}\Vert$, soit $\frac{1}{\Vert \mathbf{x}\Vert }\le \frac{\Vert A\Vert }{\Vert \mathbf{b}\Vert }$. En multipliant membre à membre ces deux inégalités entre nombres positifs,

Le conditionnement est donc le facteur d'amplification des erreurs relatives: si les données sont connues à $10^{-p}$ près et $\kappa (A)\approx 10^q$, la solution n'est garantie qu'à $10^{q-p}$ près. Avec des données en double précision ($p=16$), une matrice de conditionnement $10^{10}$ ne laisse que 6 chiffres justes. La borne est atteinte lorsque $\mathbf{b}$ est dans la direction de ${\mathbf{u}}_1$ (le plus étiré) et $\delta \mathbf{b}$ dans celle de ${\mathbf{u}}_n$ (le plus écrasé); on démontre une inégalité analogue pour une perturbation $\delta A$ de la matrice.

Conséquences pratiques

Démonstration. De $A=U\Sigma V^{\mathsf{T}}$ on tire $A^{\mathsf{T}}A=V{\Sigma }^{\mathsf{T}}{U}^{\mathsf{T}}U\Sigma V^{\mathsf{T}}=V({\Sigma }^{\mathsf{T}}\Sigma )V^{\mathsf{T}}$, avec ${\Sigma }^{\mathsf{T}}\Sigma =\mathrm{diag}({\sigma }_1^2,\dots ,{\sigma }_n^2)$. C'est une diagonalisation orthogonale de la matrice symétrique définie positive $A^{\mathsf{T}}A$, donc aussi sa SVD: ses valeurs singulières sont les ${\sigma }_i^2$ et ${\kappa }_2(A^{\mathsf{T}}A)={\sigma }_1^2/{\sigma }_n^2={\kappa }_2(A{)}^2$. $\square$

Stabilité inverse. Le conditionnement décrit la sensibilité du problème; la qualité de l'algorithme se mesure autrement. Un algorithme est dit stable au sens inverse (backward stable) si la solution calculée $\tilde{\mathbf{x}}$ est la solution exacte d'un problème voisin: $(A+\delta A)\tilde{\mathbf{x}}=\mathbf{b}$ avec $\Vert \delta A\Vert /\Vert A\Vert$ de l'ordre de $\epsilon$. La décomposition QR de Householder et la SVD sont stables au sens inverse; la décomposition LU avec pivotage partiel l'est en pratique (Wilkinson, 1961). Combinée à (12.14), la stabilité inverse donne la règle d'or: erreur relative sur la solution $\approx \kappa (A)\times \epsilon$.

Applications

Compression d'image par SVD tronquée

Une image en niveaux de gris de $m\times n$ pixels est une matrice $A$; sa SVD (12.11) l'écrit comme somme de $r$ «couches» ${\sigma }_i{\mathbf{u}}_i{\mathbf{v}}_i^{\mathsf{T}}$ d'importance décroissante. En ne gardant que les $k$ premières, on obtient la SVD tronquée

matrice de rang $k$ qui ne demande de stocker que $k(m+n+1)$ nombres au lieu de $mn$.

Nous admettons l'optimalité (la démonstration, courte mais astucieuse, se trouve dans Trefethen et Bau); les valeurs des erreurs, elles, découlent directement de la proposition 12.9 appliquée à $A-A_k={\sum }_{i>k}{\sigma }_i{\mathbf{u}}_i{\mathbf{v}}_i^{\mathsf{T}}$, dont les valeurs singulières sont ${\sigma }_{k+1},\dots ,{\sigma }_r$.

Pour une photographie de $512\times 512$ pixels, garder $k=20$ valeurs singulières demande $20\cdot (512+512+1)=20500$ nombres au lieu de $262144$, soit un taux de compression d'environ $12,8$; comme les valeurs singulières d'une image naturelle décroissent vite, le résultat est reconnaissable, quoique flou dans les détails. Le principe — représenter les données dans une base adaptée et négliger les petites composantes — est le même que dans les formats JPEG (base de cosinus) et dans la réduction de dimension des données.

Analyse en composantes principales

L'analyse en composantes principales (ACP) cherche les directions dans lesquelles un nuage de points de ${\mathbb{R}}^p$ est le plus dispersé. Soit $X\in {\mathcal{M}}_{n\times p}(\mathbb{R})$ la matrice des données centrées (une ligne par observation, une colonne par variable, la moyenne de chaque colonne ramenée à $0$). La matrice de covariance empirique est $C=\frac{1}{n-1}{X}^{\mathsf{T}}X$, et la variance des données projetées sur une direction unitaire $\mathbf{v}$ vaut ${\mathbf{v}}^{\mathsf{T}}C\mathbf{v}$. Par la SVD $X=U\Sigma V^{\mathsf{T}}$, on a $C=V\frac{{\Sigma }^{\mathsf{T}}\Sigma }{n-1}{V}^{\mathsf{T}}$: les composantes principales sont les vecteurs singuliers à droite ${\mathbf{v}}_1,{\mathbf{v}}_2,\dots$ de $X$, et la variance portée par la $i$-ième vaut ${\sigma }_i^2/(n-1)$. La preuve que ${\mathbf{v}}_1$ maximise ${\mathbf{v}}^{\mathsf{T}}C\mathbf{v}$ est celle de la proposition 12.9.

Systèmes surdéterminés en géodésie et positionnement

Le récepteur GPS d'un smartphone reçoit à chaque instant les signaux de huit à douze satellites, et cherche quatre inconnues: ses trois coordonnées et le décalage de son horloge. Les équations de distance sont non linéaires; on les linéarise autour d'une position approchée (méthode de Gauss–Newton), ce qui donne à chaque itération un système surdéterminé $A\delta \mathbf{x}\approx \mathbf{r}$ résolu par QR. La géodésie fait de même à grande échelle: la compensation par moindres carrés du réseau de triangulation suisse (cadre de référence MN95 de swisstopo) ajuste des dizaines de milliers de mesures pour quelques milliers de coordonnées inconnues. Gauss lui-même a inventé la méthode des moindres carrés pour ce type de problème, lors du levé du royaume de Hanovre.

Méthodes itératives pour les grands systèmes

Lorsque $n$ atteint le million, la décomposition LU (coût $n^3$, mémoire $n^2$) est hors de portée. Mais les matrices issues des équations aux dérivées partielles sont creuses: chaque ligne n'a que quelques coefficients non nuls, et le produit $A\mathbf{x}$ ne coûte que quelques $n$ flops. Les méthodes itératives n'utilisent que ce produit. Écrivons $A=D+L+U$ avec $D$ la diagonale (supposée sans zéro), $L$ la partie strictement inférieure et $U$ la partie strictement supérieure. Le système $A\mathbf{x}=\mathbf{b}$ équivaut à $D\mathbf{x}=\mathbf{b}-(L+U)\mathbf{x}$, ce qui suggère l'itération de Jacobi:

La méthode de Gauss–Seidel utilise immédiatement les composantes déjà mises à jour:

Nous admettons ce théorème; l'idée du critère est que l'erreur ${\mathbf{e}}^{(k)}={\mathbf{x}}^{(k)}-\mathbf{x}$ vérifie ${\mathbf{e}}^{(k+1)}=M{\mathbf{e}}^{(k)}$, donc ${\mathbf{e}}^{(k)}=M^k{\mathbf{e}}^{(0)}$, et $M^k\to 0$ si et seulement si les valeurs propres de $M$ sont de module inférieur à $1$ (chapitre 10).

Un aperçu des solveurs creux

Dans un calcul par éléments finis (chapitre 7 du cours de mécanique des structures), la matrice de rigidité d'un maillage de $10^6$ nœuds a $10^{12}$ coefficients dont seuls quelques $10^7$ sont non nuls. Deux familles de solveurs se partagent le travail. Les solveurs directs creux effectuent une décomposition LU ou de Cholesky en renumérotant les inconnues pour limiter le remplissage (les zéros de $A$ qui deviennent non nuls dans $L$); ils sont robustes et permettent de traiter de nombreux seconds membres. Les solveurs itératifs, héritiers de Jacobi et de Gauss–Seidel — méthode du gradient conjugué pour les matrices symétriques définies positives, GMRES en général — ne demandent que des produits matrice-vecteur et une mémoire proportionnelle à $n$; leur vitesse de convergence dépend directement du conditionnement, que l'on réduit par un préconditionnement (résoudre $P^{-1}A\mathbf{x}=P^{-1}\mathbf{b}$ avec $P\approx A$ facile à inverser). Le choix entre les deux, et le réglage du préconditionneur, sont le quotidien du calcul scientifique.

Synthèse

• LU: l'élimination de Gauss est la factorisation $A=LU$ ($L$ triangulaire inférieure unitaire des multiplicateurs, $U$ échelonnée); coût $2n^3/3$, puis $2n^2$ par second membre; le pivotage partiel ($PA=LU$) est indispensable en flottants; $\det A=\pm \prod u_{ii}$. Pour $A$ symétrique définie positive, Cholesky donne $A=L{L}^{\mathsf{T}}$ pour la moitié du coût, sans pivotage.
• QR: Gram–Schmidt écrit $A=QR$ ($Q$ à colonnes orthonormées, $R$ triangulaire supérieure); les réflexions de Householder $H=I-2\mathbf{v}{\mathbf{v}}^{\mathsf{T}}/{\mathbf{v}}^{\mathsf{T}}\mathbf{v}$ (symétriques, orthogonales) donnent le même résultat de façon stable. Moindres carrés: $R\hat{\mathbf{x}}=Q^{\mathsf{T}}\mathbf{b}$. L'algorithme QR ($A_{k+1}=R_k{Q}_k$) calcule les valeurs propres.
• SVD: $A=U\Sigma V^{\mathsf{T}}$ pour toute matrice; ${\sigma }_i=\sqrt{{\lambda }_i(A^{\mathsf{T}}A)}$; rotation–étirement–rotation, l'image de la sphère unité est un ellipsoïde de demi-axes ${\sigma }_i{\mathbf{u}}_i$; le rang est le nombre de ${\sigma }_i\ne 0$; $U$ et $V$ contiennent des bases orthonormées des quatre sous-espaces fondamentaux; $\Vert A{\Vert }_2={\sigma }_1$, $\Vert A{\Vert }_F^2=\sum {\sigma }_i^2$; pseudo-inverse $A^{+}=V{\Sigma }^{+}{U}^{\mathsf{T}}$ et moindres carrés de norme minimale.
• Conditionnement: $\kappa (A)=\Vert A\Vert \Vert A^{-1}\Vert$, ${\kappa }_2={\sigma }_1/{\sigma }_n$; $\Vert \delta \mathbf{x}\Vert /\Vert \mathbf{x}\Vert \le \kappa (A)\Vert \delta \mathbf{b}\Vert /\Vert \mathbf{b}\Vert$; on perd ${\log }_{10}\kappa$ chiffres; $\kappa (A^{\mathsf{T}}A)=\kappa (A{)}^2$, d'où la préférence pour QR; ne jamais inverser pour résoudre.
• Applications: SVD tronquée et théorème d'Eckart–Young (compression, débruitage), analyse en composantes principales (composantes $={\mathbf{v}}_i$, variances $={\sigma }_i^2/(n-1)$), moindres carrés en géodésie et GPS.
• Grands systèmes: méthodes itératives de Jacobi et de Gauss–Seidel (convergence si diagonale strictement dominante ou, pour Gauss–Seidel, symétrique définie positive), solveurs creux directs et itératifs préconditionnés.

Exercices

Vous pouvez afficher le corrigé directement sous chaque énoncé après avoir cherché la solution.