Raffinement itératif

Objectifs d'apprentissage

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

• Comprendre le principe du raffinement itératif
• Implémenter l'algorithme de raffinement
• Analyser quand et pourquoi cette technique fonctionne
• Identifier les situations où le raffinement est efficace
• Calculer le gain de précision attendu

Prérequis

• Conditionnement des systèmes linéaires
• Factorisation LU

Motivation : récupérer la précision perdue

Dans la leçon précédente, nous avons vu qu'un système mal conditionné peut amplifier les erreurs de manière catastrophique. Avec un conditionnement $\kappa \approx 10^k$, on perd environ $k$ chiffres significatifs.

Question : Peut-on faire mieux sans changer de méthode de résolution ?

Réponse : Oui ! Le raffinement itératif permet de récupérer une partie (voire la totalité) de la précision perdue, en réutilisant la factorisation LU déjà calculée.

Principe de l'algorithme

L'idée fondamentale

Soit $x^*$ notre solution approchée du système $Ax = b$. Définissons :

• L'erreur : $e = x - x^*$ (la vraie solution moins notre approximation)
• Le résidu : $r = b - Ax^*$ (ce qui « manque » pour satisfaire l'équation)

Ces deux quantités sont liées par une relation simple :

Pourquoi c'est utile ?

1. On ne peut pas calculer $e$ directement (on ne connaît pas $x$)
2. Mais on peut calculer $r = b - Ax^*$ (on connaît $A$, $b$ et $x^*$)
3. En résolvant $Ae = r$, on obtient une approximation de l'erreur
4. On corrige : $x^{**} = x^* + e$

L'algorithme complet

Visualisation interactive

La visualisation ci-dessous illustre le raffinement itératif sur le système mal conditionné de la leçon précédente. Cliquez sur les boutons pour voir l'évolution pas à pas avec tous les détails des calculs.

Le système utilisé

C'est le même système mal conditionné que dans la leçon sur le conditionnement :

Avec $\kappa(A) \approx 400$, on perd environ 2-3 chiffres significatifs lors de la résolution par LU.

Comment lire cette visualisation

En haut — Le graphique :

• Montre l'évolution du nombre de chiffres corrects
• Chaque point correspond à une itération

En bas — Les détails de l'itération :

Pour chaque itération, on affiche les trois quantités clés :

Ce qu'on observe

• Initial : La solution LU a environ 2 chiffres corrects (erreur $\sim 10^{-2}$)
• Iter 1 : On récupère 2 chiffres → 4 chiffres corrects
• Iter 2 : On gagne encore → 8 chiffres corrects
• Iter 3 : On atteint la précision machine (≈ 15 chiffres)

Exemple détaillé pas à pas

Prenons un système 2×2 mal conditionné pour illustrer le raffinement.

Le système

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

Conditionnement : $\kappa(A) \approx 800$ (on perd ~3 chiffres)

Étape 0 : Solution initiale (avec erreurs d'arrondi simulées)

Supposons qu'à cause des erreurs d'arrondi, notre solveur LU donne :

au lieu de $(1, 1)^T$. L'erreur initiale est :

Itération 1

Étape 1 : Calculer le résidu

Calculons $A \cdot x^{(0)}$ :

Donc :

Étape 2 : Résoudre $A \cdot e = r$

On résout :

Par élimination (ou en réutilisant LU) :

• Ligne 2 - Ligne 1 : $0{,}005 \cdot e_2 = -0{,}001015 - (-0{,}001) = -0{,}000015$
• Donc $e_2 = -0{,}003$
• De la ligne 1 : $e_1 = -0{,}001 - e_2 = -0{,}001 - (-0{,}003) = 0{,}002$

On obtient :

Étape 3 : Corriger la solution

Résultat : Après une seule itération, on a récupéré la solution exacte !

Pourquoi le raffinement fonctionne-t-il ?

Le secret : la précision du calcul du résidu

Le résidu $r = b - Ax^*$ peut être calculé avec plus de précision que la solution elle-même. Voici pourquoi :

Calcul de $x^*$ : Nécessite de résoudre un système linéaire, ce qui accumule des erreurs d'arrondi amplifiées par $\kappa(A)$.

Calcul de $r$ : C'est un simple produit matrice-vecteur suivi d'une soustraction. Les erreurs d'arrondi ne sont pas amplifiées par $\kappa$.

Analyse de convergence

Soit $e^{(k)} = x - x^{(k)}$ l'erreur à l'itération $k$. On peut montrer que :

Interprétation :

• Si $\kappa(A) \cdot \varepsilon_m < 1$, l'erreur diminue à chaque itération
• Le facteur de réduction est environ $\kappa(A) \cdot \varepsilon_m$

Exemple : Avec $\kappa = 10^8$ et $\varepsilon_m = 10^{-16}$ (double précision) :

Chaque itération réduit l'erreur d'un facteur $10^8$ ! En 2 itérations, on peut gagner 16 chiffres.

Limite de la méthode

Si $\kappa(A) > 10^{16}$, le raffinement ne converge pas — la matrice est trop mal conditionnée pour être traitée en double précision.

Amélioration : calcul du résidu en précision étendue

Pour maximiser l'efficacité du raffinement, on peut calculer le résidu en précision étendue :

Technique du « double-double »

L'idée est de calculer $r = b - Ax^*$ en accumulant les produits avec plus de précision :

import numpy as np

def residual_extended(A, b, x):
  """
  Calcule le résidu r = b - A·x avec précision étendue.
  Utilise l'algorithme de sommation compensée de Kahan.
  """
  n = len(b)
  r = np.zeros(n)

  for i in range(n):
      # Initialiser avec b[i]
      s = b[i]
      c = 0.0  # Compensation

      # Soustraire A[i,:] · x avec compensation
      for j in range(n):
          # Produit exact en deux parties
          prod = A[i, j] * x[j]

          # Sommation compensée (Kahan)
          y = -prod - c
          t = s + y
          c = (t - s) - y
          s = t

      r[i] = s

  return r

Implémentation complète

import numpy as np
from scipy.linalg import lu_factor, lu_solve

def raffinement_iteratif(A, b, max_iter=10, tol=1e-14, verbose=True):
  """
  Résout Ax = b avec raffinement itératif.

  Paramètres:
  -----------
  A : ndarray (n, n)
      Matrice du système
  b : ndarray (n,)
      Second membre
  max_iter : int
      Nombre maximum d'itérations de raffinement
  tol : float
      Tolérance sur le résidu relatif
  verbose : bool
      Afficher les informations de convergence

  Retourne:
  ---------
  x : ndarray (n,)
      Solution raffinée
  info : dict
      Informations sur la convergence
  """
  n = len(b)

  # Étape 1 : Factorisation LU (une seule fois !)
  lu, piv = lu_factor(A)

  # Étape 2 : Solution initiale
  x = lu_solve((lu, piv), b)

  # Stocker l'historique
  residuals = []

  # Étape 3 : Raffinement itératif
  for k in range(max_iter):
      # Calculer le résidu r = b - A·x
      r = b - A @ x

      # Norme relative du résidu
      res_rel = np.linalg.norm(r) / np.linalg.norm(b)
      residuals.append(res_rel)

      if verbose:
          print(f"Itération {k}: ||r||/||b|| = {res_rel:.2e}")

      # Test de convergence
      if res_rel < tol:
          if verbose:
              print(f"Convergence atteinte en {k+1} itérations")
          break

      # Résoudre A·e = r (réutilise la factorisation LU)
      e = lu_solve((lu, piv), r)

      # Corriger la solution
      x = x + e

  info = {
      'iterations': k + 1,
      'residuals': residuals,
      'final_residual': res_rel
  }

  return x, info


# ============================================================
# EXEMPLE D'UTILISATION
# ============================================================

if __name__ == "__main__":
  # Créer une matrice mal conditionnée
  n = 5

  # Matrice de Hilbert (très mal conditionnée)
  H = np.array([[1.0/(i+j+1) for j in range(n)] for i in range(n)])

  # Solution exacte connue
  x_exact = np.ones(n)

  # Second membre
  b = H @ x_exact

  # Calculer le conditionnement
  kappa = np.linalg.cond(H)
  print(f"Conditionnement κ(H) = {kappa:.2e}")
  print(f"Chiffres perdus attendus : {np.log10(kappa):.1f}")
  print()

  # Résoudre avec raffinement
  x, info = raffinement_iteratif(H, b, verbose=True)

  # Erreur finale
  erreur = np.linalg.norm(x - x_exact) / np.linalg.norm(x_exact)
  print(f"\nErreur relative finale : {erreur:.2e}")

Quand utiliser le raffinement itératif ?

Situations favorables

Avantages

• Coût faible : Chaque itération coûte $O(n^2)$ (produit matrice-vecteur + résolution triangulaire)
• Réutilise LU : La factorisation (coûteuse, $O(n^3)$) n'est faite qu'une fois
• Simple à implémenter : Quelques lignes de code supplémentaires
• Amélioration garantie : Si $\kappa \cdot \varepsilon_m < 1$, l'erreur diminue

Inconvénients

• Limité par le conditionnement : Ne fonctionne pas si $\kappa$ est trop grand
• Stockage de $A$ : On doit garder la matrice originale (pas seulement LU)
• Pas de miracle : Ne transforme pas un problème mal posé en problème bien posé

Comparaison avec et sans raffinement

Exemple numérique

Considérons la matrice de Hilbert $H_6$ avec $\kappa \approx 10^7$ :

On récupère environ 6 chiffres (ce qui correspond à $\log_{10}(\kappa) \approx 7$) !

Résumé

Le raffinement itératif est une technique simple et efficace pour améliorer la précision d'une solution obtenue par factorisation LU.

Algorithme :

1. Calculer le résidu $r = b - Ax^*$
2. Résoudre $Ae = r$ en réutilisant LU
3. Corriger $x^* \leftarrow x^* + e$
4. Répéter si nécessaire

Points clés :

• Converge si $\kappa(A) \cdot \varepsilon_m < 1$
• Chaque itération coûte $O(n^2)$ seulement
• Peut récupérer jusqu'à $\log_{10}(\kappa)$ chiffres perdus
• Le calcul du résidu en précision étendue améliore encore l'efficacité

Pour aller plus loin

La prochaine leçon présentera les méthodes itératives (Jacobi, Gauss-Seidel, SOR), une approche alternative aux méthodes directes pour les grands systèmes creux.