Conditionnement et méthodes itératives

Objectifs du chapitre

À la fin de ce chapitre, vous serez capable de:

• calculer les normes $\Vert \cdot {\Vert }_1$, $\Vert \cdot {\Vert }_2$ et $\Vert \cdot {\Vert }_{\infty }$ d'un vecteur, énoncer les inégalités d'équivalence qui les relient en dimension finie, et dire pourquoi le choix de la norme change les constantes mais pas les conclusions;
• calculer les normes matricielles induites $\Vert A{\Vert }_{\infty }$ et comme un maximum de sommes de lignes et de colonnes, et comme la plus grande valeur singulière;

Deux questions, un seul chapitre

Le chapitre 3 a fourni un algorithme qui résout $A\mathbf{x}=\mathbf{b}$ en un nombre fini d'étapes: l'élimination de Gauss et sa factorisation $LU$, avec pivot partiel. Cet algorithme est stable au sens du chapitre 1, il coûte $\frac{2}{3}{n}^3$ opérations, et il termine. Deux questions restent pourtant ouvertes, et ce sont elles qui font vivre ce chapitre.

La première est une question de fiabilité. Supposons que la factorisation se soit parfaitement passée et que le vecteur $\hat{\mathbf{x}}$ obtenu vérifie $A\hat{\mathbf{x}}=\mathbf{b}$ à $10^{-15}$ près. Que vaut $\hat{\mathbf{x}}-\mathbf{x}$? Le chapitre 1 a annoncé la réponse sans la démontrer: elle dépend d'un nombre attaché à seule, le conditionnement, et ce nombre peut valoir , ou , ou davantage. Nous allons le définir, le calculer exactement, et démontrer l'inégalité qu'il gouverne.

La seconde est une question de coût. Le facteur $n^3$ est une condamnation dès que $n$ dépasse quelques dizaines de milliers: une matrice pleine d'ordre $10^6$ demanderait $6,7\cdot 10^{17}$ opérations et $8\cdot 10^{12}$ octets de mémoire, ce qui n'existe pas. Or les grands systèmes de la pratique — discrétisation d'une équation aux dérivées partielles, réseau électrique, structure en éléments finis — sont : chaque ligne ne porte que trois, cinq ou sept coefficients non nuls. Une méthode qui ne ferait que multiplier par , sans jamais la factoriser, exploiterait cette structure. C'est l'idée des méthodes itératives, et la moitié de ce chapitre leur est consacrée.

Les deux questions se rejoignent, et c'est pourquoi elles cohabitent ici: une méthode itérative s'arrête sur un résidu, et la première moitié du chapitre explique exactement ce qu'un résidu permet de conclure.

Mesurer un vecteur: les trois normes usuelles

Pour dire qu'une erreur est petite, il faut d'abord savoir mesurer un vecteur. En dimension $1$ la valeur absolue s'impose; en dimension $n$, il y a plusieurs mesures raisonnables, et elles ne sont pas équivalentes numériquement — elles le sont seulement topologiquement, ce qui n'est pas la même chose et ce que la fin de cette section précise.

Elles ne mesurent pas la même chose et ce sont trois lectures différentes d'un même vecteur d'erreur. La norme $\Vert \cdot {\Vert }_{\infty }$ répond à «quelle est la pire composante?», et c'est celle que l'on veut quand chaque composante est une grandeur physique à garantir. La norme $\Vert \cdot {\Vert }_1$ répond à «combien d'erreur y a-t-il en tout?», et c'est celle des bilans de matière ou d'argent, où les écarts s'additionnent. La norme $\Vert \cdot {\Vert }_2$ est la longueur euclidienne, celle qui provient d'un produit scalaire, donc la seule des trois pour laquelle «orthogonal» a un sens — d'où son rôle central au chapitre 7 sur les moindres carrés.

Démonstration. Notons $m=\Vert \mathbf{x}{\Vert }_{\infty }$ et soit $k$ un indice où le maximum est atteint.

Première chaîne. On a $\Vert \mathbf{x}{\Vert }_2^2={\sum }_i{x}_i^2\ge x_k^2=m^2$, d'où ; et , d'où . Les deux bornes sont atteintes, la première par , la seconde par .

Deuxième chaîne. En développant,

puisque les termes croisés sont positifs ou nuls. Dans l'autre sens, l'inégalité de Cauchy–Schwarz appliquée à $\mathbf{x}$ et au vecteur $(1,\dots ,1)$ donne

La troisième chaîne se déduit des deux premières, ou se vérifie directement.

Cas général. Toute norme $N$ est continue pour la topologie usuelle, car $|N(\mathbf{x})-N(\mathbf{y})|\le N(\mathbf{x}-\mathbf{y})\le ({\sum }_{i}N({\mathbf{e}}_i))\Vert \mathbf{x}-\mathbf{y}{\Vert }_{\infty }$. Elle atteint donc sur la sphère unité , qui est compacte, un minimum (strictement positif par séparation) et un maximum . L'homogénéité propage ces bornes à tout .

Normes matricielles induites

Une matrice agit sur des vecteurs. La bonne façon de mesurer sa taille est donc de mesurer de combien elle peut allonger un vecteur, ce qui se dit en une formule.

L'inégalité (4.4) est la seule propriété dont nous nous servirons vraiment, et elle est la raison d'être de toute la construction. Elle se prolonge en une propriété que les normes matricielles quelconques n'ont pas toutes:

On dit qu'une norme induite est sous-multiplicative. La démonstration tient en une ligne: $\Vert AB\mathbf{x}\Vert \le \Vert A\Vert \Vert B\mathbf{x}\Vert \le \Vert A\Vert \Vert B\Vert \Vert \mathbf{x}\Vert$, puis on prend le supremum sur les $\mathbf{x}$ de norme $1$. Notons aussi que $\Vert I\Vert =1$ pour toute norme induite, ce qui n'est pas vrai de la norme de Frobenius $({\sum }_{ij}{a}_{ij}^2{)}^{1/2}$, laquelle vaut sur l'identité: la norme de Frobenius est commode à calculer mais n'est induite par aucune norme vectorielle, et nous ne l'emploierons pas.

Il reste à savoir calculer (4.3). Le miracle est que pour les deux normes extrêmes, un maximum sur une sphère se ramène à un maximum sur $n$ nombres.

Démonstration. Démontrons (4.6) en entier; (4.7) suit le même plan et (4.8) est discutée ensuite.

Posons $M={\max }_i{\sum }_j|a_{ij}|$ et soit $\mathbf{x}$ avec $\Vert \mathbf{x}{\Vert }_{\infty }=1$, c'est-à-dire pour tout . Pour toute ligne ,

En prenant le maximum sur $i$, $\Vert A\mathbf{x}{\Vert }_{\infty }\le M$, donc $\Vert A{\Vert }_{\infty }\le M$.

Pour l'égalité, il faut exhiber un vecteur qui atteint $M$. Soit $k$ une ligne réalisant le maximum. Définissons $\mathbf{z}$ par

Alors $\Vert \mathbf{z}{\Vert }_{\infty }=1$ et

d'où $\Vert A\mathbf{z}{\Vert }_{\infty }\ge M$ et donc $\Vert A{\Vert }_{\infty }\ge M$. Les deux inégalités donnent (4.6).

Pour (4.7), on écrit $A\mathbf{x}={\sum }_j{x}_j{\mathbf{a}}_j$ où ${\mathbf{a}}_j$ est la -ième colonne, d'où

et l'égalité est atteinte en prenant pour $\mathbf{x}$ le vecteur de base ${\mathbf{e}}_k$ associé à la colonne de plus grande somme. $\square$

La formule (4.8) demande davantage. Le quotient $\Vert A\mathbf{x}{\Vert }_2^2/\Vert \mathbf{x}{\Vert }_2^2$ s'écrit ${\mathbf{x}}^{\mathsf{T}}{A}^{\mathsf{T}}A\mathbf{x}/{\mathbf{x}}^{\mathsf{T}}\mathbf{x}$, c'est-à-dire le de la matrice symétrique semi-définie positive . Le théorème spectral (Algèbre linéaire) affirme que ce quotient est maximal sur le vecteur propre associé à la plus grande valeur propre, et vaut alors cette valeur propre. Les racines carrées des valeurs propres de sont les de , d'où (4.8). Nous admettons ici l'argument du quotient de Rayleigh: sa démonstration complète, et la méthode de la puissance qui le calcule, font l'objet du chapitre 5.

Une dernière relation sera indispensable dans la seconde moitié du chapitre. Le rayon spectral de $A$ est

Si $\lambda$ est une valeur propre de $A$ et $\mathbf{v}\ne 0$ un vecteur propre associé, alors, pour toute norme induite, $|\lambda |\Vert \mathbf{v}\Vert =\Vert A\mathbf{v}\Vert \le \Vert A\Vert \Vert \mathbf{v}\Vert$, d'où

Le rayon spectral est donc en dessous de toutes les normes. Il n'en est pas une lui-même: la matrice non nulle $\left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right)$ a un rayon spectral nul.

Le conditionnement d'une matrice

L'annexe A rassemble les propriétés des normes et du conditionnement dont les chapitres suivants se serviront, avec les rappels d'algèbre linéaire correspondants; nous n'en démontrons ici que ce dont ce chapitre a besoin.

Quatre propriétés élémentaires, toutes utiles:

1. $\kappa (A)\ge 1$. En effet $1=\Vert I\Vert =\Vert A{A}^{-1}\Vert \le \Vert A\Vert \Vert A^{-1}\Vert$ par (4.5).
2. $\kappa (\alpha A)=\kappa (A)$ pour tout . Multiplier une équation par ne change rien au problème, et le conditionnement le reflète: . En particulier, : la matrice a pour déterminant et pour conditionnement .

La majoration de perturbation

Voici le théorème central du chapitre. Il répond à la question: si l'on modifie légèrement le second membre — parce qu'il provient d'une mesure, d'un arrondi, ou d'un calcul antérieur — de combien la solution peut-elle bouger?

Démonstration. Les deux systèmes étant linéaires, soustrayons-les: $A\delta \mathbf{x}=\delta \mathbf{b}$, donc

par la propriété (4.4) appliquée à $A^{-1}$.

Il faut maintenant convertir cette majoration absolue en majoration relative, et c'est le second membre qui fournit le dénominateur. De $\mathbf{b}=A\mathbf{x}$ on tire, toujours par (4.4),

où l'on a utilisé $\mathbf{x}\ne 0$ (conséquence de $\mathbf{b}\ne 0$) et $\Vert \mathbf{b}\Vert >0$. En multipliant (4.13) par (4.14):

ce qui est (4.12).

Pour l'optimalité, il suffit de rendre les deux inégalités (4.13) et (4.14) simultanément exactes. Choisissons $\mathbf{x}$ réalisant le maximum de (4.3) pour $A$, c'est-à-dire tel que $\Vert A\mathbf{x}\Vert =\Vert A\Vert \Vert \mathbf{x}\Vert$, et posons $\mathbf{b}=A\mathbf{x}$: (4.14) devient une égalité. Choisissons ensuite $\delta \mathbf{b}$ réalisant le maximum de (4.3) pour $A^{-1}$, c'est-à-dire tel que $\Vert A^{-1}\delta \mathbf{b}\Vert =\Vert A^{-1}\Vert \Vert \delta \mathbf{b}\Vert$: (4.13) devient une égalité à son tour. Les deux vecteurs maximisants existent par compacité de la sphère unité, et ils sont indépendants l'un de l'autre. La borne (4.12) est donc la meilleure possible, et non une majoration confortable: .

La lecture géométrique

En dimension $2$, chaque équation $a_{i1}{x}_1+a_{i2}{x}_2=b_i$ est une droite, et la solution est le point d'intersection. Perturber la droite sans la tourner. Toute la théorie précédente se lit alors sur un dessin.

Le dessin dit exactement ce que dit le théorème 4.3, et la comparaison des deux losanges est la lecture qu'il faut emporter. Deux droites presque parallèles se coupent en un point que la moindre translation fait courir très loin le long de leur direction commune. Et un système dont les deux équations disent presque la même chose est précisément un système dont la matrice est presque singulière. Le mauvais conditionnement n'est pas un accident numérique: c'est une quasi-redondance de l'information. Quand deux capteurs mesurent presque la même combinaison de grandeurs, aucun algorithme ne pourra séparer proprement ces grandeurs, et il faut changer le dispositif expérimental, pas le solveur.

La quantification est instructive. Si les deux droites ont pour normales unitaires faisant entre elles l'angle $\theta$, on vérifie (exercice 4.2) que

avec $\theta$ en radians. Le minimum ${\kappa }_{\infty }=2$ est atteint pour $\theta =90°$: en norme infinie, le meilleur système $2\times 2$ possible est celui dont les droites sont perpendiculaires, et son conditionnement vaut $2$, pas $1$. Un angle de $1°$, soit radian, donne déjà .

Comparez cet explorateur à celui du chapitre 3, qui met en scène les deux mêmes droites mais laisse régler leurs coefficients plutôt que leur angle. Là-bas, la perturbation du second membre est quelconque, et la lecture «amplification observée» reste toujours en dessous de ${\kappa }_{\infty }$: l'inégalité (4.12) est une majoration, et une perturbation prise au hasard n'en réalise qu'une fraction. Ici, au contraire, $\delta \mathbf{b}$ est choisi dans la direction $(1,-1)$ qui sature la borne, et le rapport affiché vaut ${\kappa }_{\infty }$ exactement, quel que soit l'angle. Les deux lectures sont nécessaires: la seconde dit ce que le conditionnement peut coûter, la première ce qu'il coûte le plus souvent.

Et si c'est la matrice qui est perturbée?

Dans la pratique, $A$ aussi est connue de façon imparfaite: ses coefficients ont été mesurés, ou arrondis à l'entrée dans la machine. Le résultat correspondant, que nous admettons — sa démonstration demande le lemme de perturbation de Neumann, un exercice de série entière matricielle qui n'apporte rien de neuf ici — est le suivant: si $A$ et $A+\delta A$ sont inversibles et $\Vert A^{-1}\Vert \Vert \delta A\Vert <1$, alors

Le même facteur $\kappa (A)$ gouverne les deux sources. C'est ce qui justifie de parler du conditionnement du problème $A\mathbf{x}=\mathbf{b}$, et non d'un conditionnement par type de perturbation.

Résidu et erreur: le piège central

Nous pouvons maintenant tenir la promesse du chapitre 1.

Cette différence de statut est décisive. Le résidu se mesure, l'erreur ne se mesure pas: c'est pourquoi tout critère d'arrêt s'écrit sur le résidu, et c'est pourquoi il faut savoir ce qu'un petit résidu permet de conclure.

Démonstration. Par définition de $\mathbf{r}$, $A\hat{\mathbf{x}}=\mathbf{b}-\mathbf{r}$: le vecteur $\hat{\mathbf{x}}$ est donc la solution exacte du système dont le second membre est $\mathbf{b}-\mathbf{r}$. En posant $\delta \mathbf{b}=-\mathbf{r}$ et , on est exactement dans la situation du théorème 4.3, dont la conclusion (4.12) donne la majoration de droite.

Pour la minoration, on part de $\mathbf{r}=\mathbf{b}-A\hat{\mathbf{x}}=A\mathbf{x}-A\hat{\mathbf{x}}=-A(\hat{\mathbf{x}}-\mathbf{x})$, d'où ; et de , d'où . En divisant la première inégalité par la seconde:

ce qui est la minoration annoncée. $\square$

L'encadrement (4.18) est la formule à retenir du chapitre. Il dit deux choses symétriques, et la seconde est aussi importante que la première: si $\kappa (A)$ est petit, résidu relatif et erreur relative sont du même ordre, à un facteur $\kappa$ près des deux côtés — le résidu est alors un excellent indicateur. Si $\kappa (A)$ est grand, l'encadrement est large de ${\kappa }^2$ et ne dit plus grand-chose: un résidu minuscule est compatible avec une erreur énorme, et réciproquement.

La matrice de Hilbert, ou le mal conditionnement à l'état pur

Le mauvais conditionnement n'exige pas de coefficients pathologiques. La matrice de Hilbert d'ordre $n$,

n'a que des coefficients simples, est symétrique, définie positive, et son conditionnement est multiplié par une trentaine chaque fois que $n$ augmente de $1$. Elle apparaît naturellement: c'est la matrice de Gram de la base $1,x,x^2,\dots$ pour le produit scalaire ${\int }_0^{1}fg$, donc la matrice des équations normales de l'ajustement polynomial sur $[0,1]$ — le chapitre 7 la retrouvera, et expliquera pourquoi on n'ajuste jamais un polynôme de degré élevé par cette voie.

Aux ordres supérieurs, la situation devient caricaturale. Le tableau suivant a été calculé en arithmétique rationnelle exacte, par inversion de $H_n$ en fractions, puis arrondi pour l'affichage:

La dernière colonne est ${\log }_{10}{\kappa }_{\infty }$, c'est-à-dire l'estimation du nombre de chiffres décimaux perdus. Elle croît d'environ $1,5$ chiffre par ordre: chaque ligne ajoutée à une matrice de Hilbert coûte un chiffre et demi sur seize. À $n=11$ il n'en reste guère plus d'un, et à $n=12$, où ${\kappa }_{\infty }=4,12\cdot 10^{16}$, plus aucun.

Mettons cette prédiction à l'épreuve. Résolvons $H_n\mathbf{x}=\mathbf{b}$ en double précision par élimination de Gauss avec pivot partiel (le code du chapitre 3), avec $\mathbf{b}$ construit pour que la solution exacte soit $(1,1,\dots ,1)$:

Lisez la deuxième colonne de haut en bas: le résidu reste au niveau de la machine, invariablement, pendant que l'erreur grimpe de $10^{-14}$ à $3\cdot 10^{-4}$. L'algorithme n'a rien fait de mal — c'est même exactement le comportement d'un algorithme stable — et la troisième colonne montre que la borne (4.18) est respectée avec un facteur $5$ à $35$ de marge. À $n=10$, la solution calculée porte quatre chiffres corrects au lieu de seize. Les dernières décimales de ces trois lignes dépendent de l'ordre des opérations et ne se reproduiront pas à l'identique sur une autre implémentation; les ordres de grandeur, eux, sont robustes, et ce sont eux qui comptent.

Les méthodes itératives: décomposer pour itérer

Changeons de question. Pour $n$ grand et $A$ creuse, la factorisation $LU$ souffre d'un défaut que le chapitre 3 a nommé: le remplissage. Les zéros de $A$ ne sont pas des zéros de $L$ et de $U$, et une matrice creuse d'ordre $10^6$ peut avoir des facteurs denses, donc impossibles à stocker. Une méthode qui n'utiliserait $A$ que pour calculer des produits $A\mathbf{v}$ n'aurait pas ce problème: un produit matrice-vecteur creux coûte autant d'opérations qu'il y a de coefficients non nuls, et ne demande aucune mémoire supplémentaire.

Deux remarques valent d'être faites immédiatement.

D'abord, la solution exacte est un point fixe: si $\mathbf{x}$ vérifie $A\mathbf{x}=\mathbf{b}$, alors $M\mathbf{x}=N\mathbf{x}+\mathbf{b}$, donc $\mathbf{x}=B\mathbf{x}+\mathbf{f}$. On dit que la méthode est consistante. Cela ne garantit rien sur la convergence, exactement comme au chapitre 2 où $g(x^{\ast })=x^{\ast }$ ne garantissait pas que l'itération converge.

Ensuite, l'erreur suit une récurrence d'une simplicité totale. En soustrayant $\mathbf{x}=B\mathbf{x}+\mathbf{f}$ de (4.21) écrite avec $B$, et en posant ${\mathbf{e}}^{(k)}={\mathbf{x}}^{(k)}-\mathbf{x}$:

Toute la question de la convergence tient donc dans une seule question de théorie des matrices: quand $B^k$ tend-il vers la matrice nulle? Ce sera le théorème 4.5. Auparavant, choisissons $M$.

On écrit traditionnellement

où $D$ est la diagonale de $A$, $-E$ sa partie strictement triangulaire inférieure et $-F$ sa partie strictement triangulaire supérieure. Les deux décompositions classiques sont alors immédiates.

La méthode de Jacobi

On prend $M=D$ et $N=E+F$. Inverser $D$ est gratuit — ce sont $n$ divisions — pourvu qu'aucun coefficient diagonal ne soit nul.

La lecture de (4.24) est parlante: on résout la $i$-ième équation par rapport à sa $i$-ième inconnue, en tenant les autres pour connues. Toutes les composantes sont mises à jour à partir du seul vecteur précédent, ce qui rend la méthode trivialement parallélisable — sa seule véritable qualité aujourd'hui.

def jacobi(A, b, tol=1e-10, kmax=500):
    """Un balayage = n divisions et n(n-1) multiplications-additions, soit environ 2n^2 operations."""
    n = len(A)
    x = [0.0] * n
    nb = max(abs(v) for v in b)
    for k in range(kmax):
        r = [b[i] - sum(A[i][j] * x[j] 

Une itération coûte $n^2$ multiplications et autant d'additions pour une matrice pleine, soit environ $2n^2$ opérations, plus $2n^2$ pour le résidu si on le recalcule à chaque tour. Sur une matrice creuse à $m$ coefficients non nuls par ligne, ce coût tombe à $2mn$. À comparer aux $\frac{2}{3}{n}^3$ de la factorisation :

La méthode de Gauss–Seidel

Le calcul de $x_i^{(k+1)}$ par (4.24) utilise $x_1^{(k)},\dots ,x_{i-1}^{(k)}$ alors que viennent d'être calculés et sont, espère-t-on, meilleurs. Les utiliser tout de suite est l'idée de Gauss–Seidel.

Résoudre $M\mathbf{y}=\mathbf{c}$ avec $M$ triangulaire inférieure coûte $n^2$ opérations, mais on ne résout rien explicitement: (4.25) est la descente triangulaire, écrite ligne par ligne. Le coût par itération est donc identique à celui de Jacobi, à ceci près qu'un seul vecteur est stocké au lieu de deux, ce qui divise la mémoire par deux. En Python, la différence entre les deux méthodes tient dans une seule ligne:

def gauss_seidel(A, b, tol=1e-10, kmax=500):
    """Meme cout par iteration que Jacobi, un seul vecteur en memoire."""
    n = len(A)
    x = [0.0] * n
    nb = max(abs(v) for v in b)
    for k in range(kmax):
        r = [b[i] - sum(A[i][j] * x[j] 

L'écriture en place x[i] = ... fait que la somme de droite lit les valeurs déjà mises à jour pour $j<i$ et les anciennes pour $j>i$: c'est exactement (4.25). La contrepartie est que Gauss–Seidel dépend de l'ordre des équations, contrairement à Jacobi, et n'est pas parallélisable en l'état.

Quand l'itération converge-t-elle?

La récurrence (4.22) réduit la question à l'étude des puissances de $B$.

Démonstration. Nous démontrons $3\Rightarrow 1$ et $1\Rightarrow 2$, qui suffisent aux applications, et nous admettons $2\Rightarrow 3$.

$3\Rightarrow 1$: si $\Vert B\Vert =q<1$ pour une norme induite, la sous-multiplicativité (4.5) donne $\Vert B^k\Vert \le q^k\to 0$, donc par (4.22), et ceci quel que soit , c'est-à-dire quel que soit .

$1\Rightarrow 2$: supposons $\rho (B)\ge 1$ et soit $\lambda$ une valeur propre de module maximal, $\mathbf{v}\ne 0$ un vecteur propre associé (complexe s'il le faut). Alors $B^k\mathbf{v}={\lambda }^k\mathbf{v}$, dont la norme ne tend pas vers zéro. Donc ne tend pas vers , et l'itération initialisée en ne converge pas. Par contraposée, impose .

$2\Rightarrow 3$ est le théorème de Householder: pour tout $\epsilon >0$ il existe une norme induite telle que $\Vert B\Vert \le \rho (B)+\epsilon$. Sa démonstration passe par la forme normale de Jordan de $B$ et une mise à l'échelle des blocs; elle appartient à un cours d'analyse matricielle et n'éclairerait pas davantage le présent chapitre. Nous l'admettons.

Enfin, l'énoncé sur le taux découle de $\Vert {\mathbf{e}}^{(k)}{\Vert }^{1/k}\to \rho (B)$ (formule de Gelfand, également admise): gagner un facteur $10$ demande $k$ tel que ${\rho }^k=10^{-1}$, soit .

Le critère du théorème 4.5 est exact mais coûteux: calculer $\rho (B)$ est un problème aux valeurs propres, donc plus cher que de résoudre le système. En pratique on utilise une condition suffisante, et la plus utile de toutes est purement arithmétique.

Cette condition se lit d'un coup d'œil sur la matrice, sans aucun calcul. Elle est fréquente dans les applications: une discrétisation par différences finies d'un opérateur de diffusion avec un terme de réaction, une matrice d'admittance de réseau, une chaîne de Markov diagonalement renforcée la vérifient.

Démonstration. Commençons par observer que (4.26) impose $|a_{ii}|>{\sum }_{j\ne i}|a_{ij}|\ge 0$, donc pour tout : la matrice diagonale est inversible et l'itération de Jacobi (4.24) est bien définie.

La matrice d'itération $B_J=D^{-1}(E+F)$ a pour coefficients

Sa $i$-ième somme de ligne vaut donc, en modules,

La dominance stricte (4.26) dit exactement que $q_i<1$ pour tout $i$. Par la formule (4.6) du théorème 4.2,

le maximum d'un nombre fini de quantités strictement inférieures à $1$ étant lui-même strictement inférieur à $1$ — c'est ici, et seulement ici, que la finitude de la dimension intervient.

L'implication $3\Rightarrow 1$ du théorème 4.5 s'applique alors avec la norme $\Vert \cdot {\Vert }_{\infty }$: la récurrence de l'erreur (4.22) et la sous-multiplicativité donnent

ce qui est (4.27) et prouve la convergence pour tout ${\mathbf{x}}^{(0)}$.

Il reste l'inversibilité de $A$. Si $A\mathbf{v}=0$ avec $\mathbf{v}\ne 0$, alors $D\mathbf{v}=(E+F)\mathbf{v}$, donc $\mathbf{v}=B_J\mathbf{v}$, d'où avec : impossible pour . Le noyau de est donc réduit à zéro et est inversible.

Le même argument s'applique mot pour mot à Gauss–Seidel, avec un peu plus de travail: on montre que sous (4.26) on a aussi $\Vert B_{GS}{\Vert }_{\infty }\le q<1$, avec le même $q$ que pour Jacobi (exercice 4.6). Gauss–Seidel converge donc au moins aussi vite que ce que la borne garantit pour Jacobi. Attention toutefois: en général, ni «Jacobi converge» n'implique «Gauss–Seidel converge», ni l'inverse. Il existe des matrices où l'une converge et l'autre diverge. L'implication est vraie dans deux cas classiques: les matrices à diagonale strictement dominante, et les matrices symétriques définies positives (pour lesquelles Gauss–Seidel converge toujours, ce qui n'est pas le cas de Jacobi).

Les disques de la figure 4.2 sont ceux du théorème de Gershgorin, démontré au chapitre 5: toute valeur propre d'une matrice $B$ appartient à la réunion des disques de centre $b_{ii}$ et de rayon ${\sum }_{j\ne i}|b_{ij}|$. Appliqué à , dont la diagonale est nulle, il redonne exactement le théorème 4.6 — les disques sont centrés à l'origine et de rayons , donc . C'est le même calcul que (4.27), vu autrement:

La relaxation: accélérer en dosant le pas

Gauss–Seidel calcule, pour chaque composante, une valeur $x_i^{GS}$ censée améliorer $x_i^{(k)}$. L'idée de la relaxation est de ne pas s'arrêter à cette valeur, mais d'aller plus ou moins loin dans la direction de la correction.

L'écriture (4.28) se lit comme une moyenne pondérée entre l'ancienne valeur et la valeur de Gauss–Seidel: $x_i^{(k+1)}=(1-\omega )x_i^{(k)}+\omega x_i^{GS}$. Avec on délibérément la correction proposée, dans l'idée que Gauss–Seidel corrige toujours un peu trop peu. Deux résultats encadrent le choix de .

Démonstration de (a). Le déterminant d'un produit est le produit des déterminants, et celui d'une matrice triangulaire le produit de sa diagonale. Or $\frac{1}{\omega }D-E$ est triangulaire inférieure de diagonale $\frac{1}{\omega }D$, et $(\frac{1}{\omega }-1)D+F$ est triangulaire supérieure de diagonale . Donc

Le déterminant étant le produit des $n$ valeurs propres, ${\prod }_i|{\lambda }_i|=|1-\omega {|}^n$, et le plus grand module est au moins la moyenne géométrique: $\rho (B_{\omega })\ge |1-\omega |$. Si , c'est-à-dire si ou , alors et le théorème 4.5 interdit la convergence.

Les points (b) et (c) sont admis: le premier repose sur une identité d'énergie ${\mathbf{v}}^{\mathsf{T}}A\mathbf{v}$ un peu longue, le second sur l'analyse de la relation entre les valeurs propres de $B_{\omega }$ et celles de $B_J$ pour les matrices cohéremment ordonnées (relation de Young), qui relève d'un cours spécialisé. Ce qu'il faut retenir de (4.30) est sa conséquence pratique, spectaculaire.

Le système de la figure 4.3 est $T\mathbf{x}=\mathbf{b}$ avec $T$ tridiagonale d'ordre $n=10$, de diagonale $2$ et de sous-diagonales $-1$ — c'est la matrice de la dérivée seconde discrétisée, la matrice modèle de toute la discipline. Ses rayons spectraux se calculent exactement:

et (4.30) donne

Le tableau des itérations effectivement comptées, pour un résidu relatif inférieur à $10^{-8}$, confirme tout:

Trois observations, dont la troisième est celle qui sert.

D'abord, le gain est considérable: $38$ itérations contre $197$, pour un coût par itération identique et une ligne de code de différence. Contre Jacobi, qui demande $391$ balayages sur le même système, le facteur est de $10$.

Ensuite, la courbe est très pointue au voisinage de l'optimum: de $\omega =1,5$ à $\omega =1,5604$, soit un écart de $4\%$, le nombre d'itérations passe de $56$ à $38$. Choisir $\omega$ au hasard perd l'essentiel du bénéfice.

Enfin, et c'est la règle de terrain, la courbe n'est pas symétrique: à distance égale de l'optimum, $\omega =1,5$ coûte $56$ itérations et $\omega =1,62$ en coûterait environ $43$. Pour $\omega >{\omega }_{\text{opt}}$ on a exactement $\rho (B_{\omega })=\omega -1$, une droite de pente ; en dessous, décroît beaucoup plus vite. , jamais l'inverse.

Directe ou itérative? Et quand s'arrêter?

Le comptage qui décide

Comparons honnêtement, en opérations, sur une matrice pleine d'ordre $n$:

• factorisation $LU$ puis deux substitutions: $\frac{2}{3}{n}^3+2n^2$ opérations, une fois pour toutes;
• $K$ itérations de Jacobi ou de Gauss–Seidel: environ $2K{n}^2$ opérations.

L'itératif l'emporte si $2K{n}^2<\frac{2}{3}{n}^3$, c'est-à-dire si $K<n/3$. Pour $n=1000$, il faudrait converger en moins de itérations — ce qui arrive — mais on aurait quand même stocké la matrice pleine, et l'on n'aurait rien gagné en mémoire.

Tout change sur une matrice creuse. Soit $m$ le nombre moyen de coefficients non nuls par ligne ($m=3$ pour une tridiagonale, $m=5$ pour le laplacien 2D à cinq points, $m=7$ en 3D). Alors:

• une itération coûte $2mn$ opérations et aucune mémoire supplémentaire;
• la factorisation, elle, subit le remplissage, et son coût dépend de la structure. Pour la matrice du laplacien sur une grille $N\times N$, donc $n=N^2$ inconnues et une largeur de bande $N$, une factorisation de Cholesky en bande coûte environ $\frac{1}{2}n{N}^2=\frac{1}{2}{n}^2$ opérations et demande cases mémoire.

Chiffrons, avec une tolérance de $10^{-8}$ et $\omega$ optimal (pour lequel le nombre d'itérations croît comme $N$):

Le croisement se produit vers $N\approx 60$, et l'écart s'ouvre ensuite: le coût direct croît comme $n^2$, le coût itératif comme $n^{3/2}$. À un million d'inconnues, SOR est $17$ fois moins cher en opérations et $200$ fois moins cher en mémoire — et c'est la mémoire qui décide en premier, car un milliard de cases, c'est $8$ gigaoctets rien que pour les facteurs.

Ce sont des comptages d'opérations, pas des temps mesurés. Deux implémentations du même comptage peuvent différer d'un facteur $10$ en secondes selon l'accès mémoire, et ce cours ne prétend pas prédire des secondes.

Les critères d'arrêt

Une méthode itérative doit décider seule quand s'arrêter, et le seul objet dont elle dispose est le résidu. Le critère standard est donc

Trois remarques, dans l'ordre d'importance.

(1) Il est relatif, et il doit l'être. Un critère absolu $\Vert {\mathbf{r}}^{(k)}\Vert \le \text{tol}$ dépend des unités: multiplier $\mathbf{b}$ par $1000$ — passer de mètres à millimètres — changerait le nombre d'itérations. Le rapport (4.31) est sans dimension.

(2) Il borne l'erreur à travers $\kappa$, et pas autrement. Le théorème 4.4 donne, quand (4.31) est satisfait,

Choisir $\text{tol}=10^{-8}$ sur un système de conditionnement $10^6$ ne garantit que deux chiffres. La tolérance se choisit en fonction du conditionnement, ou au moins de son ordre de grandeur estimé.

(3) Le critère sur l'incrément est un piège. Il est tentant d'arrêter sur $\Vert {\mathbf{x}}^{(k+1)}-{\mathbf{x}}^{(k)}\Vert \le \text{tol}$, qui ne coûte rien. Mais l'incrément et l'erreur sont reliés par

Si $\rho (B)$ est proche de $1$ — et c'est exactement le cas où la méthode est lente, donc où l'on est tenté d'arrêter tôt — le facteur $\Vert (I-B{)}^{-1}\Vert$ est énorme, et un incrément minuscule accompagne une erreur considérable. Pour $\rho =0,999$, l'erreur peut valoir mille fois l'incrément. Un critère sur $\Vert {\mathbf{x}}^{(k+1)}-{\mathbf{x}}^{(k)}\Vert$ n'est pas une borne sur l'erreur; c'est le même avertissement qu'au chapitre 2 pour les itérations scalaires, et il s'aggrave ici parce que est souvent très proche de .

En production, on combine les trois: un critère relatif sur le résidu, un plafond d'itérations, et un drapeau de non-convergence — c'est exactement ce que construit l'exercice 4.8.

Méthodes de gradient pour les matrices symétriques définies positives

Jacobi, Gauss–Seidel et la relaxation traitent toutes les matrices de la même façon, et la meilleure d'entre elles demande un paramètre $\omega$ que l'on ne connaît pas. Or une grande partie des systèmes de la pratique ont une structure de plus: la matrice est symétrique définie positive (SDP, définition A.6 de l'annexe A). C'est le cas de la matrice $T$ de la figure 4.3, des matrices de rigidité des éléments finis, des matrices de diffusion, et de la matrice $A^{\mathsf{T}}A$ des équations normales du chapitre 7. Pour elles, résoudre le système revient à minimiser une énergie, et cette reformulation donne une famille de méthodes entièrement différente, sans paramètre à régler, dont le gradient conjugué est le représentant le plus utilisé du calcul scientifique.

Résoudre, c'est minimiser

Démonstration. Posons $\mathbf{e}=\mathbf{x}-{\mathbf{x}}^{\ast }$ et développons $\phi ({\mathbf{x}}^{\ast }+\mathbf{e})$ en utilisant la symétrie, qui donne ${\mathbf{e}}^{\mathsf{T}}A{\mathbf{x}}^{\ast }={\mathbf{x}}^{\ast \mathsf{T}}A\mathbf{e}$:

puisque $A{\mathbf{x}}^{\ast }=\mathbf{b}$. Comme $\Vert \mathbf{e}{\Vert }_A>0$ pour $\mathbf{e}\ne 0$, $\phi (\mathbf{x})>\phi ({\mathbf{x}}^{\ast })$ pour tout . Le gradient se lit sur le même développement écrit autour d'un quelconque: , dont la partie linéaire en est (Analyse II, chapitre 4).

Géométriquement, les courbes de niveau de $\phi$ sont des ellipses (des ellipsoïdes en dimension $n$) centrées en ${\mathbf{x}}^{\ast }$, d'axes dirigés selon les vecteurs propres de $A$ et de demi-axes proportionnels à $1/\sqrt{{\lambda }_i}$. Leur — le rapport du grand au petit axe — vaut donc , puisque pour une matrice SDP (propriété 4 qui suit la définition 4.3). Le conditionnement, que la première moitié du chapitre voyait comme une sensibilité aux perturbations, devient ici une : un grand, c'est une vallée longue et étroite, et c'est là que les méthodes de descente peinent.

La plus profonde descente

La façon la plus naturelle de descendre une vallée est de suivre la pente. En ${\mathbf{x}}^{(k)}$, la direction de plus forte décroissance de $\phi$ est $-\nabla \phi ({\mathbf{x}}^{(k)})={\mathbf{r}}^{(k)}$: le résidu lui-même. On avance donc dans la direction du résidu, d'un pas ${\alpha }_k$ choisi pour minimiser le long de cette droite. La fonction est un trinôme du second degré, minimal en

La dernière relation se vérifie en appliquant $A$ à la deuxième; elle évite de recalculer $\mathbf{b}-A{\mathbf{x}}^{(k+1)}$, de sorte qu'une itération ne coûte qu'un seul produit matrice-vecteur, $A{\mathbf{r}}^{(k)}$, plus deux produits scalaires et deux mises à jour de vecteurs: $2mn+8n$ opérations sur une matrice creuse à $m$ coefficients par ligne.

def plus_profonde_descente(A, b, tol=1e-8, kmax=10000):
    n = len(b)
    x = [0.0] * n
    r = b[:]                                      # r = b - A x avec x = 0
    nb = max(abs(v) for v in b)
    for k in range(kmax):
        if max(abs(v) 

Le choix optimal de ${\alpha }_k$ a une conséquence qu'il faut voir: la dérivée du trinôme s'annule en ${\alpha }_k$, ce qui s'écrit ${\mathbf{r}}^{(k+1)\mathsf{T}}{\mathbf{r}}^{(k)}=0$. : chaque pas s'arrête là où la droite de descente est tangente à une courbe de niveau, et le pas suivant repart à angle droit. Dans une vallée étroite, cela produit le zigzag de la figure 4.4.

Démonstration. L'idée est de comparer le pas optimal à un pas fixe bien choisi. Pour un $\alpha$ quelconque, le point ${\mathbf{x}}^{(k)}+\alpha {\mathbf{r}}^{(k)}$ a pour erreur ${\mathbf{e}}^{(k)}+\alpha {\mathbf{r}}^{(k)}=(I-\alpha A){\mathbf{e}}^{(k)}$, puisque . Par (4.33), minimiser sur la droite revient à minimiser la norme de l'énergie de l'erreur; le pas optimal fait donc que n'importe quel pas fixe:

Développons ${\mathbf{e}}^{(k)}={\sum }_i{c}_i{\mathbf{v}}_i$ sur une base orthonormée de vecteurs propres de $A$ (théorème spectral, annexe A). Alors $\Vert {\mathbf{e}}^{(k)}{\Vert }_A^2={\sum }_i{\lambda }_i{c}_i^2$ et , d'où

Choisissons $\alpha =2/({\lambda }_{\min }+{\lambda }_{\max })$. La fonction $\lambda \mapsto |1-\alpha \lambda |$ est maximale sur $[{\lambda }_{\min },{\lambda }_{\max }]$ à l'une des deux extrémités, où elle vaut dans les deux cas .

Le taux $(\kappa -1)/(\kappa +1)=1-2/(\kappa +1)$ est proche de $1$ dès que $\kappa$ est grand: gagner un chiffre demande de l'ordre de $\kappa$ itérations. La plus profonde descente est lente exactement quand le système est mal conditionné, et le zigzag en est la raison géométrique: dans une vallée étroite, la pente pointe presque perpendiculairement à la direction où il faudrait aller.

Le gradient conjugué

Le zigzag vient de ce que chaque pas défait en partie le travail des précédents: après avoir minimisé dans la direction ${\mathbf{r}}^{(0)}$, puis dans ${\mathbf{r}}^{(1)}$, on n'est plus au minimum dans la direction ${\mathbf{r}}^{(0)}$, et il faudra y revenir. L'idée du gradient conjugué, due à Hestenes et Stiefel (1952), est de choisir des directions qui ne se gâchent pas mutuellement.

Si l'on dispose de $n$ directions ${\mathbf{p}}_0,\dots ,{\mathbf{p}}_{n-1}$ deux à deux $A$-conjuguées, elles forment une base, orthogonale pour $⟨\cdot ,\cdot {⟩}_A$, et l'erreur initiale s'y décompose: . Minimiser le long de détermine exactement le coefficient — c'est tout l'intérêt de l'orthogonalité —, si bien qu'après minimisations unidimensionnelles on a la solution exacte. Le gradient conjugué construit ces directions au fur et à mesure, chacune à partir du résidu courant corrigé d'un multiple de la direction précédente:

Le coefficient ${\beta }_k$ est exactement celui qui rend ${\mathbf{p}}^{(k+1)}$ conjuguée à ${\mathbf{p}}^{(k)}$. Le miracle — c'est le mot qu'emploient les manuels — est qu'elle l'est alors automatiquement à toutes les directions précédentes, sans qu'on ait à les stocker.

def gradient_conjugue(A, b, tol=1e-8, kmax=1000):
    n = len(b)
    x = [0.0] * n
    r = b[:]
    p = r[:]
    rr = sum(v * v for v in r)
    nb = max(abs(v) for v in

Une itération coûte un produit $A{\mathbf{p}}^{(k)}$, deux produits scalaires et trois mises à jour de vecteurs, soit $2mn+10n$ opérations sur une matrice creuse — à peine plus que la plus profonde descente, et la mémoire se limite à quatre vecteurs. La matrice n'intervient que par des produits matrice-vecteur, comme dans toutes les méthodes itératives de ce chapitre.

Démonstration. Les points 1 et 2 s'établissent ensemble par récurrence sur $k$, en utilisant à chaque étape la définition de ${\alpha }_k$ (qui rend ${\mathbf{r}}^{(k+1)}$ orthogonal à ${\mathbf{r}}^{(k)}$), celle de ${\beta }_k$ (qui rend conjuguée à ) et la symétrie de pour les indices plus anciens. Le calcul est élémentaire mais long, et nous l'admettons; on le trouve dans Quarteroni, Sacco et Saleri, ou dans Saad (voir les références). Le point 3 en découle: par récurrence, et appartiennent à , et l'erreur est -orthogonale à , qui engendrent : c'est la caractérisation de la meilleure approximation pour le produit scalaire de l'énergie.

La terminaison finie, elle, se démontre en une ligne à partir du point 1. Si ${\mathbf{r}}^{(0)},\dots ,{\mathbf{r}}^{(n-1)}$ sont tous non nuls, ce sont $n$ vecteurs non nuls deux à deux orthogonaux de ${\mathbb{R}}^n$, donc une base orthogonale; ${\mathbf{r}}^{(n)}$, orthogonal à tous, est nul. Donc .

Le point 3 dit plus que la terminaison: il dit que le gradient conjugué fait, à chaque itération, le mieux possible avec l'information qu'il a recueillie. Il en découle une borne qui fait apparaître la racine carrée du conditionnement.

Nous admettons ce théorème, dont voici l'idée. Tout élément de ${\mathbf{x}}^{(0)}+{\mathcal{K}}_k$ a une erreur de la forme $P(A){\mathbf{e}}^{(0)}$, où $P$ est un polynôme de degré au plus $k$ avec . Par le point 3 et le même développement spectral que dans la démonstration du théorème 4.9, pour tel . La plus profonde descente revenait à prendre à chaque pas; le gradient conjugué choisit implicitement le meilleur polynôme de degré , et le meilleur polynôme pour ce problème de minimax est un polynôme de déplacé — les mêmes polynômes dont le chapitre 6 tire ses nœuds d'interpolation. Leur croissance hors de donne (4.37).

Comparez (4.35) et (4.37): $(\kappa -1)/(\kappa +1)$ contre $(\sqrt{\kappa }-1)/(\sqrt{\kappa }+1)$. Pour , les taux valent et : gagner six chiffres demande environ itérations dans le premier cas et dans le second, en ignorant le facteur 2 de (4.37). , et c'est ce qui en fait la méthode de référence pour les grands systèmes SDP.

Le préconditionnement. La borne (4.37) dit que tout se joue sur $\kappa$. Si l'on dispose d'une matrice $M$ symétrique définie positive, proche de $A$ et facile à inverser, on applique le gradient conjugué au système équivalent $M^{-1}A\mathbf{x}=M^{-1}\mathbf{b}$, réécrit de façon symétrique: c'est le gradient conjugué préconditionné. Une itération demande en plus la résolution d'un système $M\mathbf{z}=\mathbf{r}$, et la vitesse est gouvernée par le conditionnement de au lieu de celui de . Le choix de est un arbitrage du même type que celui de la décomposition des méthodes de Jacobi et de Gauss–Seidel. Le plus simple est la diagonale, — qui ne change rien sur , dont la diagonale est constante, et c'est un bon rappel qu'un préconditionneur se juge sur le problème. Les plus utilisés sont la , qui calcule les facteurs du chapitre 3 en refusant tout remplissage, et les méthodes multigrilles, pour lesquelles le nombre d'itérations ne dépend plus de sur les matrices de discrétisation. Le principe est toujours le même: .

Synthèse

• Les trois normes $\Vert \cdot {\Vert }_1$, $\Vert \cdot {\Vert }_2$ et $\Vert \cdot {\Vert }_{\infty }$ sont équivalentes en dimension finie, avec des constantes en $\sqrt{n}$ ou (théorème 4.1): les conclusions qualitatives ne dépendent pas du choix, les chiffres si. , y compris sur le conditionnement.

Exercices

Vous pouvez afficher le corrigé directement sous chaque énoncé après avoir cherché la solution.

Références

• Quarteroni, A., Sacco, R. et Saleri, F., Méthodes numériques — algorithmes, analyse et applications, Springer, Milan, chap. 3 (normes et conditionnement) et chap. 4 (méthodes itératives pour les systèmes linéaires).
• Rappaz, J. et Picasso, M., Introduction à l'analyse numérique, Presses polytechniques et universitaires romandes, Lausanne, chap. 3.
• Trefethen, L. N. et Bau, D., Numerical Linear Algebra, SIAM, Philadelphie, leçons 3 à 4 (normes matricielles) et 12 (conditionnement d'un système linéaire).
• Higham, N. J., Accuracy and Stability of Numerical Algorithms, 2e éd., SIAM, Philadelphie, chap. 7 (perturbation des systèmes linéaires) et chap. 9.
• Saad, Y., Iterative Methods for Sparse Linear Systems, 2e éd., SIAM, Philadelphie, chap. 4 (méthodes de relaxation et leur convergence), chap. 5 et 6 (plus profonde descente, espaces de Krylov et gradient conjugué) et chap. 9 et 10 (préconditionnement).
• Hestenes, M. R. et Stiefel, E., «Methods of conjugate gradients for solving linear systems», Journal of Research of the National Bureau of Standards, vol. 49, 1952 (l'article fondateur du gradient conjugué).
• Burden, R. L. et Faires, J. D., Numerical Analysis, Cengage, Boston, chap. 7.