2. Préconditionnement

Objectifs

• Ajouter un préconditionneur
• Comparer les performances

Principe

Le nombre d’itérations nécessaire pour la méthode du Gradient Conjugué peut être diminué. Nous savons que la convergence de la méthode dépendant du conditionnement de la matrice du système linéaire. Or, pour la matrice du Laplacien, ce conditionnement tend vers l’infini quand la taille $N$ du système tend également vers l’infini. Cela nous amène à étudier la notion de préconditionnement.

Prenons un système linéaire $Ax = b$. Préconditionner ce système revient à multiplier le système par une matrice $P^{-1}$ et à résoudre le système linéaire suivant: $$ P^{-1}Ax = P^{-1}b. $$ Le vecteur $x$ est également solution du système original. La matrice $P$ doit être choisie de telle sorte que :

1. Elle soit “facilement” inversible, au sens où le coût de calcul de $P^{-1}$ doit rester faible
2. Le conditionnement de la matrice $P^{-1}A$ soit le plus proche possible de 1 ($P^{-1}A$ doit être proche de l’identité).
3. Si l’on souhaite utiliser la méthode du Gradient Conjugué pour résoudre le nouveau système, $P^{-1}A$ soit toujours symétrique et définie positive.

Notez que le meilleur préconditionneur possible, au sens du critère 2, est bien entendu $P^{-1} = A^{-1}$, mais cela ne respecterait pas le critère 1.

Préconditionneur simple : Jacobi

Un préconditionneur très simple est celui de Jacobi, $P_J := \text{diag}(A)$. Dans le cas d’une matrice à coefficients diagonaux identiques, comme la matrice du Laplacien, ce préconditionneur n’a aucune influence sur la convergence.

Matrice test

Nous proposons de travailler sur une autre matrice symétrique définie positive $H_n$, de taille $n\times n$, telle que les seuls coefficients non nuls sont donnés par : $$ \forall i = 0,\ldots, n-1, \qquad \left\{ \begin{array}{r c l l} H_n(i,i) &=& i+1,\\\
H_n(i,i+2) &=& 1 & \text{ si } i+2 \leq n\\\
H_n(i,i-2) &=& 1 & \text{ si } i-2 \geq 0 \end{array} \right. $$

Gradient Conjugué Précondtionné

Pseudo-code

Le gradient conjugué préconditionné est modifié ainsi ((z,r) représente le produit scalaire en z et r)

x_0 = 0
r_0 = b - A*x_0
z_0 = P^{-1}*r0
p_0 = z_0
k = 0
while (|r| / |b| > tolerance && k < nmax)
  k++
  alpha_k = (r_k, zk)/(p_k, A*p_k)
  x_{k+1} = x_k + alpha_k * p_k
  r_{k+1} = r_k - alpha_k*(A*p_k)
  z_{k+1} = P^{-1}*r_{k+1}
  beta_k  = (z_{k+1}, r_{k+1})/(z_k,r_k)
  p_{k+1} = z_{k+1} + beta_k*p_k
end
return x_{k+1}

Historique de convergence

Soit $b_n$ le vecteur de taille $n$ rempli de 1, nous considérons le problème linéaire suivant \begin{equation} \label{eq:hnsystem} H_n x = b_n, \end{equation} et sa version préconditionnée par Jacobi (où $P_n = \text{diag}(H_n)$) \begin{equation} \label{eq:pnhnsystem} P_n^{-1}H_n x = P_n^{-1}b_n. \end{equation}