D´ecompositions LU et Choleski

Ce document présente les décompositions LU et de Choleski, deux méthodes fondamentales en algèbre linéaire pour la résolution de systèmes linéaires et l'inversion de matrices. Destiné aux étudiants préparant l’épreuve d’informatique de l’École Polytechnique, il expose les définitions, conditions d’existence, algorithmes de calcul, ainsi que des exemples et démonstrations, en utilisant le langage Maple.

La décomposition LU

Soit A une matrice carrée inversible d’ordre n sur un corps K. Une décomposition LU de A est une écriture :

A = LU

où L est une matrice triangulaire inférieure à diagonale unité (tous les coefficients diagonaux valent 1) et U est une matrice triangulaire supérieure inversible (ses coefficients diagonaux sont non nuls).

On note a_ij, l_ij et u_ij les coefficients généraux de A, L et U respectivement. Alors :

a_ij = ∑_k=1^min(i,j) l_ik u_kj

En particulier :

• Pour i ≤ j : u_ij = a_ij − ∑_k=1ⁱ⁻¹ l_ik u_kj
• Pour i > j : l_ij = (1 / u_jj) (a_ij − ∑_k=1^j−1 l_ik u_kj)

Ces relations permettent de calculer progressivement les coefficients de L et U, colonne par colonne.

Condition d'existence et unicité

La décomposition LU existe et est unique si et seulement si tous les mineurs principaux de A sont non nuls. Les mineurs principaux ∆_k sont les déterminants des sous-matrices carrées extraites des k premières lignes et colonnes de A.

Procédure LU1 : calcul direct de L et U

On peut traduire les formules précédentes en une procédure algorithmique (ici en pseudo-code Maple) :

LU1 := proc(A::array, L::name, U::name, n::integer)
local i, j, k, c;
L := matrix(n, n, 0);
U := matrix(n, n, 0);
for j to n do
  for i to j do
    c := A[i, j];
    for k to i-1 do
      c := c - L[i, k] * U[k, j];
    od;
    U[i, j] := c;
  od;
  L[j, j] := 1;
  for i from j+1 to n do
    c := A[i, j];
    for k to j-1 do
      c := c - L[i, k] * U[k, j];
    od;
    L[i, j] := c / U[j, j];
  od;
od;
end:

Cette procédure calcule explicitement L et U à partir de A.

Procédure LU2 : décomposition "in situ"

Pour économiser la mémoire, on peut stocker simultanément L et U dans la même matrice B, définie par :

• b_ij = u_ij si i ≤ j
• b_ij = l_ij si i > j

La formule de calcul devient :

Pour j de 1 à n, et pour i de 1 à n :

b_ij = a_ij − ∑_k=1^min(i,j)-1 b_ik b_kj

Si i > j, alors :

b_ij ← b_ij / b_jj

Cette méthode permet de modifier la matrice A directement en B, évitant la création de matrices supplémentaires.

Procédure LU2 (pseudo-code)

LU2 := proc(A::array, n)
local i, j, k, p;
for j to n do
  for i to n do
    for k to min(i, j) - 1 do
      A[i, j] := A[i, j] - A[i, k] * A[k, j];
    od;
    if i > j then
      A[i, j] := A[i, j] / A[j, j];
    fi;
  od;
od;
end:

Décomposition avec permutation des lignes : P A = LU

La décomposition LU classique nécessite que les pivots (coefficients diagonaux de U) soient non nuls. Si un pivot est nul ou trop petit, on effectue des échanges de lignes. Cela revient à décomposer une matrice P A, où P est une matrice de permutation.

Une matrice de permutation P correspond à une permutation σ de {1, ..., n} telle que :

p_ij = δ_{σ(i), j}

où δ est le symbole de Kronecker.

On a la propriété suivante :

Proposition : Pour toute matrice inversible A, il existe une matrice de permutation P telle que P A admet une décomposition LU.

Procédure LU avec pivot partiel

Cette procédure recherche à chaque étape le pivot de valeur absolue maximale dans la colonne courante, échange les lignes correspondantes, puis effectue la décomposition in situ.

LU := proc(A::array, s::name, n::integer)
local i, j, k, p, q, t;
s := array(1..n);
for i to n do s[i] := i od;
for j to n do
  p := A[j, j];
  i := j;
  for k from j+1 to n do
    if abs(A[k, j]) > abs(p) then
      p := A[k, j];
      i := k;
    fi;
  od;
  if p = 0 then
    ERROR("matrice non inversible");
  fi;
  if i <> j then
    for k to n do
      t := A[j, k];
      A[j, k] := A[i, k];
      A[i, k] := t;
    od;
    t := s[j];
    s[j] := s[i];
    s[i] := t;
  fi;
  for k from j+1 to n do
    A[k, j] := A[k, j] / p;
  od;
  for i from j+1 to n do
    for k from j+1 to n do
      A[i, k] := A[i, k] - A[i, j] * A[j, k];
    od;
  od;
od;
end:

Résolution d’un système linéaire AX = B à l’aide de la décomposition LU

Soit P la matrice de permutation associée à la décomposition P A = LU. Le système AX = B est équivalent à :

LU X = P B

On pose Y = U X, alors :

• Résoudre LY = P B par substitution avant (système triangulaire inférieur)
• Résoudre U X = Y par substitution arrière (système triangulaire supérieur)

Les formules sont :

• Pour i de 1 à n :

y_i = b_σ(i) − ∑_j=1ⁱ⁻¹ l_ij y_j

• Pour i de n à 1 :

x_i = (y_i − ∑_j=i+1ⁿ u_ij x_j) / u_ii

Procédure syst (pseudo-code)

syst := proc(A::array, B::array, s::array, n::integer)
local X, i, j, t;
X := array(1..n);
for i to n do
  X[i] := B[s[i]];
od;
for i from 2 to n do
  t := 0;
  for j to i-1 do
    t := t + A[i, j] * X[j];
  od;
  X[i] := X[i] - t;
od;
for i from n to 1 by -1 do
  t := 0;
  for j from i+1 to n do
    t := t + A[i, j] * X[j];
  od;
  X[i] := (X[i] - t) / A[i, i];
od;
RETURN(X);
end:

Exemple numérique

Considérons la matrice :

2	−3	1	−1
−2	2	−3	2
4	−9	−2	3
−2	5	5	−4

et le vecteur second membre B = (1, 6, 5, 2). La décomposition LU avec permutation permet de résoudre le système AX = B avec une solution X = (−17, −10, −2, −7).

On remarque que la méthode LU avec pivot partiel est plus précise que la méthode LU2 sans permutation, notamment lorsque certains coefficients sont très petits.

Inversion de matrice par décomposition LU

Pour inverser une matrice A, on résout n systèmes linéaires AX = B où B est chaque colonne de la matrice identité.

inv := proc(A::array, n::integer)
local B, M, AA, s, i, j;
M := array(1..n, 1..n);
AA := copy(A);
LU(AA, 's', n);
for j to n do
  B := array(1..n, 0);
  B[j] := 1;
  B := syst(AA, B, s, n);
  for i to n do
    M[i, j] := B[i];
  od;
od;
RETURN(eval(M));
end:

Cette méthode est efficace et donne un résultat identique à l’inversion intégrée de Maple.

La décomposition de Choleski

Considérons une matrice S réelle, carrée d’ordre n, symétrique et définie positive. Cela signifie que :

• Toutes les valeurs propres de S sont strictement positives.
• Pour tout vecteur non nul X, tX S X > 0.
• Il existe une matrice inversible M telle que S = M tM.

La décomposition de Choleski affirme qu’il existe une unique matrice triangulaire inférieure L à coefficients diagonaux strictement positifs telle que :

S = L tL

Calcul des coefficients de L

En identifiant les coefficients, on obtient :

• Pour i = j :

l_jj² = s_jj − ∑_k=1^j−1 l_jk²

• Pour i > j :

l_ij = (1 / l_jj) (s_ij − ∑_k=1^j−1 l_ik l_jk)

Ces formules permettent de calculer L colonne par colonne.

Procédure choleski (pseudo-code)

choleski := proc(S::array, n::integer)
local L, i, j, k, t;
L := matrix(n, n, 0);
for j to n do
  t := S[j, j];
  for k to j-1 do
    t := t - L[j, k]^2;
  od;
  if t <= 0 then
    ERROR("matrice non définie positive");
  fi;
  L[j, j] := sqrt(t);
  for i from j+1 to n do
    t := S[i, j];
    for k to j-1 do
      t := t - L[i, k] * L[j, k];
    od;
    L[i, j] := t / L[j, j];
  od;
od;
RETURN(eval(L));
end:

Exemple

Soit la matrice symétrique :

1	−2	4
−2	13	−11
4	−11	21

La décomposition de Choleski donne :

1	0	0
−2	3	0
4	−1	2

et on vérifie que L tL = S.

Glossaire des termes clés

• Décomposition LU : Écriture d’une matrice carrée A sous la forme A = LU, où L est triangulaire inférieure à diagonale unité et U triangulaire supérieure inversible.
• Mineurs principaux : Déterminants des sous-matrices carrées extraites des k premières lignes et colonnes de A.
• Matrice de permutation : Matrice obtenue en permutant les lignes d’une matrice identité selon une permutation σ.
• Pivot : Coefficient diagonal utilisé pour la division lors de la décomposition LU.
• Décomposition in situ : Méthode de calcul où L et U sont stockées dans la même matrice, remplaçant progressivement A.
• Décomposition de Choleski : Écriture d’une matrice symétrique définie positive S sous la forme S = L tL, avec L triangulaire inférieure à diagonale strictement positive.
• Symétrie : Propriété d’une matrice S telle que S = tS.
• Définie positive : Matrice S telle que tX S X > 0 pour tout vecteur non nul X.
• Substitution avant : Méthode de résolution d’un système triangulaire inférieur LY = B.
• Substitution arrière : Méthode de résolution d’un système triangulaire supérieur U X = Y.

Points clés à retenir

• La décomposition LU existe et est unique si tous les mineurs principaux de la matrice sont non nuls.
• La décomposition LU peut être réalisée in situ, stockant L et U dans la même matrice.
• La décomposition LU avec permutation des lignes (pivot partiel) garantit l’existence de la décomposition même si certains pivots sont nuls ou petits.
• La résolution d’un système linéaire AX = B par LU se fait en deux étapes : résolution de LY = P B puis U X = Y.
• La décomposition de Choleski est spécifique aux matrices symétriques définies positives et permet une factorisation unique S = L tL.
• Les algorithmes présentés sont adaptés à une implémentation dans Maple avec des procédures claires et efficaces.