Factorisation LU — Résolution et applications

Objectifs d'apprentissage

À la fin de cette leçon, vous serez en mesure de :

Résoudre un système linéaire en utilisant une factorisation LU existante
Comprendre et appliquer le stockage compact de L et U
Implémenter l'algorithme complet de factorisation LU avec résolution
Calculer le déterminant d'une matrice à partir de sa factorisation LU
Intégrer le pivotage dans la factorisation LU

Prérequis

Factorisation LU — Principes
Substitution avant et arrière

Résolution en deux étapes

Principe

Une fois la factorisation $A = L \cdot U$ obtenue, résoudre $A \cdot x = b$ devient :

L \cdot U \cdot x = b

On pose $y = U \cdot x$ et on résout en deux étapes :

🚨

Méthode de résolution

Substitution avant : Résoudre $L \cdot y = b$ pour trouver $y$
Substitution arrière : Résoudre $U \cdot x = y$ pour trouver $x$

Chaque étape coûte $O(n^2/2)$ , soit un total de $O(n^2)$ par système résolu.

Exemple détaillé

Résolvons le système $A \cdot x = b$ avec :

A = \begin{pmatrix} 2 & 1 & 1 \\ 4 & 3 & 3 \\ 8 & 7 & 9 \end{pmatrix}, \quad b = \begin{pmatrix} 4 \\ 10 \\ 24 \end{pmatrix}

La factorisation LU (calculée dans la leçon précédente) est :

L = \begin{pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 4 & 3 & 1 \end{pmatrix}, \quad U = \begin{pmatrix} 2 & 1 & 1 \\ 0 & 1 & 1 \\ 0 & 0 & 2 \end{pmatrix}

Étape 1 : Résoudre $L \cdot y = b$

\begin{pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 4 & 3 & 1 \end{pmatrix} \begin{pmatrix} y_1 \\ y_2 \\ y_3 \end{pmatrix} = \begin{pmatrix} 4 \\ 10 \\ 24 \end{pmatrix}

Par substitution avant :

$y_1 = 4$
$y_2 = 10 - 2 \cdot 4 = 2$
$y_3 = 24 - 4 \cdot 4 - 3 \cdot 2 = 24 - 16 - 6 = 2$

Donc $y = (4, 2, 2)^T$ .

Étape 2 : Résoudre $U \cdot x = y$

\begin{pmatrix} 2 & 1 & 1 \\ 0 & 1 & 1 \\ 0 & 0 & 2 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 4 \\ 2 \\ 2 \end{pmatrix}

Par substitution arrière :

$x_3 = 2/2 = 1$
$x_2 = (2 - 1 \cdot 1)/1 = 1$
$x_1 = (4 - 1 \cdot 1 - 1 \cdot 1)/2 = 2/2 = 1$

Solution : $x = (1, 1, 1)^T$

Avantage pour plusieurs seconds membres

Quand on doit résoudre plusieurs systèmes avec la même matrice :

Méthode	1 système	k systèmes
Gauss répété	$n^3/3$	$k \cdot n^3/3$
LU + résolution	$n^3/3 + n^2$	$n^3/3 + k \cdot n^2$

Pour $n = 1000$ et $k = 100$ systèmes :

Gauss répété : ≈ 33 milliards d'opérations
LU : ≈ 433 millions d'opérations (76× plus rapide !)

Stockage compact

Principe

Puisque $L$ a des 1 sur la diagonale (qu'on n'a pas besoin de stocker) et $U$ a des 0 sous la diagonale, on peut stocker les deux matrices dans un seul tableau de la taille de $A$ :

\begin{pmatrix} u_{11} & u_{12} & u_{13} \\ l_{21} & u_{22} & u_{23} \\ l_{31} & l_{32} & u_{33} \end{pmatrix}

Algorithme avec stockage compact

lu_compact.pypython

import numpy as np

def lu_decomposition_compact(A):
  """
  Factorisation LU avec stockage compact (in-place).
  La matrice A est modifiée pour contenir L (sous la diagonale)
  et U (diagonale et au-dessus).
  """
  n = len(A)
  LU = A.astype(float).copy()

  for k in range(n - 1):
      # Vérifier le pivot
      if abs(LU[k, k]) < 1e-12:
          raise ValueError(f"Pivot nul à la position ({k}, {k})")

      # Calculer les multiplicateurs et les stocker dans L
      for i in range(k + 1, n):
          LU[i, k] = LU[i, k] / LU[k, k]

          # Mettre à jour U
          for j in range(k + 1, n):
              LU[i, j] -= LU[i, k] * LU[k, j]

  return LU

def lu_solve_compact(LU, b):
  """
  Résout Ax = b en utilisant la factorisation LU compacte.
  """
  n = len(b)
  x = b.astype(float).copy()

  # Substitution avant (Ly = b)
  for i in range(1, n):
      for j in range(i):
          x[i] -= LU[i, j] * x[j]

  # Substitution arrière (Ux = y)
  for i in range(n - 1, -1, -1):
      for j in range(i + 1, n):
          x[i] -= LU[i, j] * x[j]
      x[i] /= LU[i, i]

  return x

# Exemple
A = np.array([[2, 1, 1],
            [4, 3, 3],
            [8, 7, 9]], dtype=float)
b = np.array([4, 10, 24], dtype=float)

LU = lu_decomposition_compact(A)
print("LU compact =")
print(LU)

x = lu_solve_compact(LU, b)
print(f"\nSolution : x = {x}")

Factorisation LU avec pivotage

Nécessité du pivotage

Comme pour l'élimination de Gauss, le pivotage est nécessaire pour :

Éviter les divisions par zéro (pivot nul)
Améliorer la stabilité numérique (pivot petit)

Factorisation PA = LU

Avec pivotage partiel, on obtient :

P \cdot A = L \cdot U

où $P$ est une matrice de permutation (enregistrant les échanges de lignes).

lu_pivotage.pypython

import numpy as np

def lu_decomposition_pivot(A):
  """
  Factorisation LU avec pivotage partiel.
  Retourne LU (compact), et perm (vecteur de permutation).
  """
  n = len(A)
  LU = A.astype(float).copy()
  perm = list(range(n))  # Permutation initiale [0, 1, 2, ...]

  for k in range(n - 1):
      # Pivotage partiel : trouver le plus grand pivot
      max_idx = k + np.argmax(np.abs(LU[k:, k]))

      if max_idx != k:
          # Échanger les lignes
          LU[[k, max_idx]] = LU[[max_idx, k]]
          perm[k], perm[max_idx] = perm[max_idx], perm[k]

      if abs(LU[k, k]) < 1e-12:
          raise ValueError("Matrice singulière")

      # Élimination
      for i in range(k + 1, n):
          LU[i, k] /= LU[k, k]
          for j in range(k + 1, n):
              LU[i, j] -= LU[i, k] * LU[k, j]

  return LU, perm

def lu_solve_pivot(LU, perm, b):
  """
  Résout Ax = b avec la factorisation PA = LU.
  """
  n = len(b)

  # Appliquer la permutation à b
  pb = np.array([b[perm[i]] for i in range(n)], dtype=float)

  # Substitution avant
  for i in range(1, n):
      for j in range(i):
          pb[i] -= LU[i, j] * pb[j]

  # Substitution arrière
  for i in range(n - 1, -1, -1):
      for j in range(i + 1, n):
          pb[i] -= LU[i, j] * pb[j]
      pb[i] /= LU[i, i]

  return pb

# Exemple avec pivot nécessaire
A = np.array([[0, 1, 2],
            [1, 2, 3],
            [2, 1, 1]], dtype=float)
b = np.array([5, 10, 6], dtype=float)

LU, perm = lu_decomposition_pivot(A)
print(f"Permutation : {perm}")
x = lu_solve_pivot(LU, perm, b)
print(f"Solution : x = {x}")

Calcul du déterminant

Le déterminant de $A$ se calcule facilement à partir de la factorisation LU :

\det(A) = \det(L) \cdot \det(U) = 1 \cdot \prod_{i=1}^{n} u_{ii} = \prod_{i=1}^{n} u_{ii}

Avec pivotage, il faut tenir compte du signe :

\det(A) = (-1)^s \prod_{i=1}^{n} u_{ii}

où :

$n$ est la dimension de la matrice
$s$ est le nombre d'échanges de lignes effectués lors du pivotage
$u_{ii}$ sont les éléments diagonaux de $U$

determinant_lu.pypython

def determinant_lu(LU, num_swaps=0):
  """
  Calcule le déterminant à partir de la factorisation LU.
  """
  diag_product = 1.0
  for i in range(len(LU)):
      diag_product *= LU[i, i]
  return ((-1) ** num_swaps) * diag_product

Avantages et inconvénients de LU

Avantages	Inconvénients
Efficace pour plusieurs seconds membres	N'exploite pas la symétrie de la matrice
Stockage compact possible	Pivotage plus complexe à gérer
Calcul du déterminant facile	Pas adapté aux matrices creuses
Préserve la structure de la matrice	Requiert mineurs principaux non nuls (sans pivotage)

Résumé

La factorisation LU permet de résoudre efficacement plusieurs systèmes linéaires :

Factorisation : $A = L \cdot U$ (une seule fois, $O(n^3/3)$ )
Résolution : $L \cdot y = b$ puis $U \cdot x = y$ (pour chaque $b$ , $O(n^2)$ )

Points clés :

Stockage compact : L et U dans une seule matrice
Pivotage : $P \cdot A = L \cdot U$ pour la stabilité
Déterminant : produit de la diagonale de U

Pour aller plus loin

La prochaine leçon abordera les systèmes rectangulaires (plus d'équations que d'inconnues ou l'inverse) et la méthode des moindres carrés pour trouver la meilleure solution approximative.

Factorisation LU — Résolution et applications

Objectifs d'apprentissage

Prérequis

Résolution en deux étapes

Principe

Méthode de résolution

Exemple détaillé

Étape 1 : Résoudre L⋅y=bL \cdot y = bL⋅y=b

Étape 2 : Résoudre U⋅x=yU \cdot x = yU⋅x=y

Avantage pour plusieurs seconds membres

Stockage compact

Principe

Algorithme avec stockage compact

Factorisation LU avec pivotage

Nécessité du pivotage

Factorisation PA = LU

Calcul du déterminant

Avantages et inconvénients de LU

Résumé

Pour aller plus loin

Étape 1 : Résoudre $L \cdot y = b$

Étape 2 : Résoudre $U \cdot x = y$