Annexe A — Factorisation de Cholesky

Introduction

La factorisation de Cholesky est une variante de la factorisation LU, spécialement conçue pour les matrices symétriques définies positives. Elle exploite la structure particulière de ces matrices pour offrir une méthode plus efficace et numériquement stable.

Matrices symétriques définies positives

Définition

Une matrice est symétrique définie positive (SDP) si :

1. est symétrique :
2. Pour tout vecteur non nul :

Propriétés importantes

Les matrices SDP possèdent plusieurs propriétés remarquables :

Exemples de matrices SDP

Les matrices SDP apparaissent naturellement dans de nombreux contextes :

• Matrices de covariance en statistiques
• Matrices de rigidité en mécanique des structures
• Matrices dans les moindres carrés
• Matrices laplaciennes en analyse numérique

Principe de la factorisation

L'idée fondamentale

Pour une matrice SDP , on cherche une matrice triangulaire inférieure telle que :

Contrairement à la factorisation LU classique où avec , ici on exploite la symétrie : la matrice triangulaire supérieure est simplement la transposée de .

Structure de la décomposition

Pour une matrice 3×3 :

Formules de calcul

Dérivation des formules

En développant le produit , on obtient :

Pour les éléments diagonaux () :

D'où :

Pour les éléments hors diagonale () :

D'où :

Algorithme

import numpy as np

def cholesky(A):
  """
  Factorisation de Cholesky : A = L·L^T

  Paramètres:
      A : matrice symétrique définie positive (n x n)

  Retourne:
      L : matrice triangulaire inférieure
  """
  n = len(A)
  L = np.zeros((n, n))

  for i in range(n):
      # Élément diagonal
      somme = sum(L[i, k]**2 for k in range(i))
      val = A[i, i] - somme

      if val <= 0:
          raise ValueError(f"Matrice non définie positive (étape {i})")

      L[i, i] = np.sqrt(val)

      # Éléments sous la diagonale
      for j in range(i + 1, n):
          somme = sum(L[i, k] * L[j, k] for k in range(i))
          L[j, i] = (A[j, i] - somme) / L[i, i]

  return L

# Exemple
A = np.array([[4, 2, 2],
            [2, 5, 1],
            [2, 1, 6]], dtype=float)

L = cholesky(A)
print("L =")
print(L)
print("\nVérification L·L^T =")
print(L @ L.T)

Exemple détaillé

Factorisons la matrice :

Étape 1 : Première colonne

Étape 2 : Deuxième colonne

Étape 3 : Troisième colonne

Résultat

Vérification

Résolution de systèmes

Pour résoudre avec :

Étape 1 : Résoudre (substitution avant)

Étape 2 : Résoudre (substitution arrière)

def resoudre_cholesky(A, b):
  """
  Résout A·x = b où A est symétrique définie positive.
  """
  L = cholesky(A)
  n = len(b)

  # Substitution avant : L·y = b
  y = np.zeros(n)
  for i in range(n):
      y[i] = (b[i] - sum(L[i, k] * y[k] for k in range(i))) / L[i, i]

  # Substitution arrière : L^T·x = y
  x = np.zeros(n)
  for i in range(n - 1, -1, -1):
      x[i] = (y[i] - sum(L[j, i] * x[j] for j in range(i + 1, n))) / L[i, i]

  return x

Avantages et inconvénients

Avantages

Inconvénients

Comparaison avec LU

Variante : Cholesky LDL^T

Pour éviter les racines carrées, on peut utiliser la factorisation :

où est triangulaire inférieure avec des 1 sur la diagonale, et est diagonale.

Cette variante est parfois préférée car elle évite les racines carrées, au prix d'une matrice supplémentaire.

Applications pratiques

En NumPy/SciPy

import numpy as np
from scipy.linalg import cholesky, cho_solve, cho_factor

# Matrice symétrique définie positive
A = np.array([[4, 2, 2],
            [2, 5, 1],
            [2, 1, 6]], dtype=float)
b = np.array([1, 2, 3], dtype=float)

# Factorisation
L = cholesky(A, lower=True)
print("L =\n", L)

# Résolution efficace
c, low = cho_factor(A)
x = cho_solve((c, low), b)
print("Solution x =", x)

Cas d'usage typiques

1. Moindres carrés : est toujours SDP
2. Krigeage/Interpolation : Matrices de covariance
3. Optimisation : Matrices hessiennes (si définies positives)
4. Éléments finis : Matrices de rigidité

Résumé

La factorisation de Cholesky est la méthode de choix pour les matrices symétriques définies positives :

• Forme :
• Condition : symétrique et définie positive
• Complexité : (2× plus rapide que LU)
• Stabilité : Toujours stable, pas de pivotage nécessaire
• Détection : Échec de la racine carrée → matrice non SDP