Décomposition LU

La méthode du pivot de Gauss n'est pas exactement programmée comme elle a été présentée. Il y a plusieurs raisons à cela, dont la principale est le problème de la précision numérique.

Voici un système de deux équations à deux inconnues, dépendant du paramètre $\varepsilon \neq 0$.

$$\left\{\begin{array}{rcl} \varepsilon x+y&=&1\ [2ex] x+y&=&2 \en... ...] y&=&(2-\frac{1}{\varepsilon })/(1-\frac{1}{\varepsilon }) \end{array}\right.$$

Voici le même système, après avoir échangé les deux équations.

$$\left\{\begin{array}{rcl} x+y&=&2\ [2ex] \varepsilon x+y&=&1 \en... ...varepsilon )\ [2ex] y&=&(1-2\varepsilon )/(1-\varepsilon ) \end{array}\right.$$

Les deux solutions sont évidemment les mêmes. Pourtant, si $\varepsilon$ est très petit en valeur absolue, les deux calculs ne sont pas du tout équivalents numériquement : diviser par un petit nombre, ou multiplier par un grand nombre, augmente les erreurs d'approximation. Telles que nous les avons présentées, les permutations de lignes et de colonnes servent à assurer que les pivots restent non nuls. La plupart des systèmes que l'on rencontre en pratique ont une solution unique : ce sont des systèmes de $n$ équations à $n$ inconnues, de rang $n$. En général, on peut leur appliquer la méthode du pivot de Gauss sans rencontrer de pivot nul. Mais on utilise quand même les permutations de lignes et de colonnes, pour faire en sorte qu'à chaque étape, le pivot soit le plus grand possible en valeur absolue.

Permuter les lignes d'une matrice, revient à la multiplier à gauche par une matrice de permutation. Une matrice de permutation est la matrice de passage de la base $(b_1,\ldots,b_n)$ à la base $(b_{\sigma(1)},\ldots,b_{\sigma(n)})$, où $\sigma$ est une bijection de $\{1,\ldots,n\}$ dans lui-même. Ses coefficients d'ordre $(\sigma(i),i)$ valent 1, les autres 0. Permuter les colonnes d'une matrice, revient à la multiplier à droite par une autre matrice de permutation. En permutant les lignes et les colonnes, on remplace la matrice $A$ par la matrice $P_1AP_2$ où $P_1$ et $P_2$ sont deux matrices de permutation.

Dans sa version la plus courante, l'algorithme ne considère que des permutations de lignes : il remplace donc la matrice $A$ par $PA$, où $P$ est une matrice de permutation. Une fois choisi l'ordre dans lequel on traite les lignes, la $i$-ième étape de la méthode consiste à ajouter aux lignes d'indice $i+1,i+2,\ldots,n$ la $i$-ième ligne multipliée par un certain coefficient. Cela revient à multiplier à gauche par une matrice du type suivant.

$$\left(\begin{array}{cccccc} 1&0&\cdots&\cdots&\cdots&0\\ 0&\ddot... ...dots&&\vdots&&\ddots&0\\ 0&\cdots&\lambda_{m,i}&0&\cdots&1 \end{array}\right)$$

Le produit de ces matrices, pour $i$ allant de $1$ à $m$ est la matrice ci-dessous.

$$\left(\begin{array}{cccccc} 1&0&\cdots&\cdots&\cdots&0\\ \lambda... ...\lambda_{m,1}&\cdots&\lambda_{m,i}&\cdots&\lambda_{m,m-1}&1 \end{array}\right)$$

Son inverse est encore une matrice du même type : trianglaire inférieure avec des $1$ sur la diagonale. On la note $L$ (pour «lower triangular»). Le produit $L^{-1}PA$ est une matrice triangulaire supérieure, que l'on note $U$ pour «upper triangular» : $U$ est la forme échelonnée de $A$.

$$L^{-1}PA=U\Longleftrightarrow PA=LU\;.$$

La décomposition LU de la matrice $A$ est la donnée des trois matrices $P,L,U$ telles que $PA=LU$. Si on doit résoudre le système $Ax=b$, on le transformera en deux systèmes triangulaires, un de matrice $L$, l'autre de matrice $U$.

$$Ax=b \Longleftrightarrow PAx=Pb \Longleftrightarrow LUx=Pb \Longleftrightarrow \left\{\begin{array}{l} Ly=Pb\ [2ex] Ux=y \end{array}\right.$$

Il arrive fréquemment que l'on ait à résoudre successivement de nombreux systèmes linaires ayant tous la même matrice $A$, mais des seconds membres différents. Calculer au préalable la décomposition LU de $A$ réduit de beaucoup le temps de calcul. Pour certaines matrices qui reviennent souvent dans les calculs, la décomposition LU figure dans les bibliothèques de codes, et elle est chargée en mémoire avant le début du calcul.