3. Factorisation LU

Objectifs

1. Calculer la factorisation LU d’une matrice
2. Résoudre le système linéaire une fois la factorisation effectuée

Principe

Cette méthode permet de transformer une matrice carré $A$ en un produit d’une matrice triangulaire inférieur $L$ et d’une matrice triangulaire supérieur $U$. Cette décomposition permet notamment de résoudre des problèmes d’algèbre linéaire du type \begin{equation} \label{eq:Axb} AX=b \iff LUX = B \end{equation} où $X$ et $B$ sont respectivement les vecteurs solution et second membre. Au final, la solution $X$ est obtenu par deux résolutions successives : $$ X = U^{-1}(L^{-1}B). $$

Ainsi, comme $L$ et $U$ sont triangulaires respectivement inférieur et supérieur, les trois étapes pour résoudre le système \eqref{eq:Axb} sont :

1. Calculer la factorisation LU
2. Résoudre un système linéaire triangulaire inférieur (avec des 1 sur la diagonale)
3. Réoudre un système linéaire triangulaire supérieur

C’est parfait, nous avons déjà implémenté les fonctions de résolution, il ne nous manque plus que le calcul de la factorisation LU !

Factorisations partielle et complète

Factorisation partielle

Notons $a_{i,j}$ le coefficient $(i,j)$ de la matrice $A$. Nous allons tout d’abord faire une factorisation partielle de la matrice

\begin{equation} A=\begin{pmatrix} a_{0,0} & A_{0,1} \\\
A_{1,0} & A_{1,1} \end{pmatrix} = \begin{pmatrix} 1 & 0 \\\
L_{1,0} & I \end{pmatrix} \begin{pmatrix} u_{0,0} & U_{0,1} \\\
0 & S_{1,1} \end{pmatrix} \label{eq:factorisation_partielle} \end{equation} où $I$ est la matrice identité, les $A_{I,J}$ sont des sous-blocs de $A$ (notons que $A_{0,0} = a_{0,0}$ est un coefficient). Le bloc $S_{1,1}=A_{1,1}-A_{1,0}A_{0,0}^{-1}A_{0,1}$ est appelé le complément de Schur.

Factorisation complète

Le lien entre factorisation partielle et factorisation complète est donné par le théorème suivant :

Ce théorème nous dit que dès lors qu’on arrive à décomposer un bloc de la diagonale $A_{0,0}$ sous forme $LU$, nous n’avons plus qu’à calculer $L_{1,0}$, $U_{0,1}$ et $S_{1,1}$ puis on cherche la décomposition $LU$ de $S_{1,1}$. Autrement dit, si nous disposons d’une fonction permettant de réaliser une factorisation partielle d’une matrice donnée, nous pouvons envisager un algorithme itératif pour obtenir la factorisation complète de la matrice.

Algorithme

Principe

Pour obtenir la factorisation complète, un algorithme itératif possible consiste à appliquer la factorisation partiellement successivement sur les compléments de Schur $S_{k,k}$ : $$ A = L^{(0)} U^{(0)}= \ldots = L^{(k)} U^{(k)} = \ldots = L^{(N-1)} U^{(N-1)}. $$ où les matrices $L^{(k)}$ et $U^{(k)}$ sont obtenues à la $k^{\text{eme}}$ itération. La petite animation suivante montre la forme de ces matrices dans le cas d’une taille N=5 :

Pseudo code

L = 0;
U = 0;
S = A;
for k =0:N-1
  // Pivot
  pivot = S(k,k)
  // Colonne de L
  L(k,k) = 1;
  for i = k+1:N-1
    L(i,k) = S(i,k) / pivot;
  // Ligne de U
  U(k,k) = S(k,k);
  for j = k+1:N-1
    U(k,j) = S(k,j);
  // Complément de Schur
  for i = k+1:N-1
    for j = k+1:N-1
      S(i,j) = S(i,j) - L(i,k)*U(k,j);

Factorisation sur place

Plutôt que de stocker 3 matrices L, U et S, dont on sait que cela coûte très cher, on remarque que l’on peut se passer de …:

• … la matrice S en modifiant directement U : le bloc $U_{k,k}$ (en “bas à droite”) contiendra le complément de Schur
• … la matrice L en la stockant dans U et en supprimant son terme diagonal (qui vaut 1 et peut donc devenir “implicite”)
• … la matrice U et travailler directement dans A

Cela donne le pseudo-code suivant :

S = A;
L = 0;
U = 0;
for k =0:N-1
  // Pivot
  pivot = S(k,k)
  // Colonne de L
  L(k,k) = 1;
  for i = k+1:N-1
    L(i,k) = S(i,k) / pivot;
  // Ligne de U
  U(k,k) = S(k,k);
  for j = k+1:N-1
    U(k,j) = S(k,j);
  // Complément de Schur
  for i = k+1:N-1
    for j = k+1:N-1
      S(i,j) = S(i,j) - L(i,k)*U(k,j);

L = 0;
U = A;
for k =0:N-1
  // Pivot
  pivot = U(k,k)
  // Colonne de L
  L(k,k) = 1;
  for i = k+1:N-1
    L(i,k) = U(i,k) / pivot;
  // Ligne de U
  U(k,k) = U(k,k);
  for j = k+1:N-1
    U(k,j) = U(k,j);
  // Complément de Schur
  for i = k+1:N-1
    for j = k+1:N-1
      U(i,j) = U(i,j) - L(i,k)*U(k,j);

U = A;
for k =0:N-1
// Pivot
pivot = U(k,k)
// Colonne de L
// U(k,k) = 1;
for i = k+1:N-1
U(i,k) = U(i,k) / pivot;
// Ligne de U
U(k,k) = U(k,k);
for j = k+1:N-1
U(k,j) = U(k,j);
// Complément de Schur
for i = k+1:N-1
for j = k+1:N-1
U(i,j) = U(i,j) - U(i,k)*U(k,j);

for k =0:N-1
// Pivot
pivot = A(k,k)
// Colonne de L
// A(k,k) = 1;
for i = k+1:N-1
A(i,k) = A(i,k) / pivot;
// Ligne de U
A(k,k) = A(k,k);
for j = k+1:N-1
A(k,j) = A(k,j);
// Complément de Schur
for i = k+1:N-1
for j = k+1:N-1
A(i,j) -= A(i,k)*A(k,j);

Implémentation en C++

void Matrice::decomp_LU();

Validation