Chapitre 4
Normes vectorielles et matricielles, nombre de conditionnement, résidu et erreur, méthodes de Jacobi, de Gauss-Seidel et de relaxation, critères de convergence.
À la fin de ce chapitre, vous serez capable de:
Le chapitre 3 a fourni un algorithme qui résout en un nombre fini d'étapes: l'élimination de Gauss et sa factorisation , avec pivot partiel. Cet algorithme est stable au sens du chapitre 1, il coûte 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 obtenu vérifie à près. Que vaut ? 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 est une condamnation dès que dépasse quelques dizaines de milliers: une matrice pleine d'ordre demanderait opérations et 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 creux: 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.
Pour dire qu'une erreur est petite, il faut d'abord savoir mesurer un vecteur. En dimension la valeur absolue s'impose; en dimension , 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 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 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 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 et soit un indice où le maximum est atteint.
Première chaîne. On a , 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 à et au vecteur donne
La troisième chaîne se déduit des deux premières, ou se vérifie directement.
Cas général. Toute norme est continue pour la topologie usuelle, car . 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 .
Un calcul en dimension renvoie un vecteur d'erreur dont vous savez seulement que . Que pouvez-vous garantir sur la plus grande composante de , c'est-à-dire sur ?
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: , puis on prend le supremum sur les de norme . Notons aussi que pour toute norme induite, ce qui n'est pas vrai de la norme de Frobenius , 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 nombres.
Démonstration. Démontrons (4.6) en entier; (4.7) suit le même plan et (4.8) est discutée ensuite.
Posons et soit avec , c'est-à-dire pour tout . Pour toute ligne ,
En prenant le maximum sur , , donc .
Pour l'égalité, il faut exhiber un vecteur qui atteint . Soit une ligne réalisant le maximum. Définissons par
Alors et
d'où et donc . Les deux inégalités donnent (4.6).
Pour (4.7), on écrit où est la -ième colonne, d'où
et l'égalité est atteinte en prenant pour le vecteur de base associé à la colonne de plus grande somme.
La formule (4.8) demande davantage. Le quotient s'écrit , c'est-à-dire le quotient de Rayleigh 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 valeurs singulières 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 est
Si est une valeur propre de et un vecteur propre associé, alors, pour toute norme induite, , 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 a un rayon spectral nul.
Écrivez norme_inf_vecteur, norme_inf_matrice et conditionnement_inf. Le programme affiche à quatre décimales, puis et pour le système témoin du cours, à une décimale. Vous devez retrouver et .
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:
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: , donc
par la propriété (4.4) appliquée à .
Il faut maintenant convertir cette majoration absolue en majoration relative, et c'est le second membre qui fournit le dénominateur. De on tire, toujours par (4.4),
où l'on a utilisé (conséquence de ) et . 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 réalisant le maximum de (4.3) pour , c'est-à-dire tel que , et posons : (4.14) devient une égalité. Choisissons ensuite réalisant le maximum de (4.3) pour , c'est-à-dire tel que : (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: le facteur est réellement atteint par certaines données.
En dimension , chaque équation est une droite, et la solution est le point d'intersection. Perturber translate 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 , on vérifie (exercice 4.2) que
avec en radians. Le minimum est atteint pour : en norme infinie, le meilleur système possible est celui dont les droites sont perpendiculaires, et son conditionnement vaut , pas . Un angle de , soit radian, donne déjà .
Les deux droites se coupent toujours au point (1, 1), et le second membre subit toujours la même perturbation relative. Resserrez l'angle: les droites deviennent presque parallèles, le conditionnement grimpe comme 2/θ, et le point d'intersection perturbé, le point plein au bout du segment, s'éloigne d'autant. Surveillez l'avant-dernière ligne des résultats: le rapport entre l'erreur relative et la perturbation relative ne bouge pas de κ∞, parce que cette perturbation-là sature exactement la borne du théorème 4.3.
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 : l'inégalité (4.12) est une majoration, et une perturbation prise au hasard n'en réalise qu'une fraction. Ici, au contraire, est choisi dans la direction qui sature la borne, et le rapport affiché vaut 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.
Deux droites du plan se coupent sous un angle de . En utilisant l'approximation de (4.15), avec exprimé en radians, estimez le conditionnement en norme infinie du système correspondant.
Dans la pratique, 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 et sont inversibles et , alors
Le même facteur gouverne les deux sources. C'est ce qui justifie de parler du conditionnement du problème , et non d'un conditionnement par type de perturbation.
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 , : le vecteur est donc la solution exacte du système dont le second membre est . En posant 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 , d'où ; et de , d'où . En divisant la première inégalité par la seconde:
ce qui est la minoration annoncée.
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 est petit, résidu relatif et erreur relative sont du même ordre, à un facteur près des deux côtés — le résidu est alors un excellent indicateur. Si est grand, l'encadrement est large de et ne dit plus grand-chose: un résidu minuscule est compatible avec une erreur énorme, et réciproquement.
Le mauvais conditionnement n'exige pas de coefficients pathologiques. La matrice de Hilbert d'ordre ,
n'a que des coefficients simples, est symétrique, définie positive, et son conditionnement est multiplié par une trentaine chaque fois que augmente de . Elle apparaît naturellement: c'est la matrice de Gram de la base pour le produit scalaire , donc la matrice des équations normales de l'ajustement polynomial sur — 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 en fractions, puis arrondi pour l'affichage:
| chiffres perdus | ||||
|---|---|---|---|---|
| 3 | 11/6 | 408 | 748 | 2,9 |
| 4 | 25/12 | 13 620 | 28 375 | 4,5 |
| 5 | 137/60 | 413 280 | 943 656 | 6,0 |
| 6 | 49/20 | 11 865 420 | 29 070 279 | 7,5 |
| 8 | 761/280 | 10,5 | ||
| 10 | 7381/2520 |
La dernière colonne est , c'est-à-dire l'estimation du nombre de chiffres décimaux perdus. Elle croît d'environ chiffre par ordre: chaque ligne ajoutée à une matrice de Hilbert coûte un chiffre et demi sur seize. À il n'en reste guère plus d'un, et à , où , plus aucun.
Mettons cette prédiction à l'épreuve. Résolvons en double précision par élimination de Gauss avec pivot partiel (le code du chapitre 3), avec construit pour que la solution exacte soit :
| erreur relative observée | résidu relatif observé | borne résidu | |
|---|---|---|---|
| 3 | |||
| 6 | |||
| 10 |
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 à . 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 à de marge. À , 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.
Vous obtenez une solution approchée d'un système linéaire et vous voulez savoir combien de chiffres lui faire confiance. Remettez les étapes dans l'ordre.
Glissez les éléments pour les mettre dans le bon ordre
Changeons de question. Pour grand et creuse, la factorisation souffre d'un défaut que le chapitre 3 a nommé: le remplissage. Les zéros de ne sont pas des zéros de et de , et une matrice creuse d'ordre peut avoir des facteurs denses, donc impossibles à stocker. Une méthode qui n'utiliserait que pour calculer des produits 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 vérifie , alors , donc . On dit que la méthode est consistante. Cela ne garantit rien sur la convergence, exactement comme au chapitre 2 où ne garantissait pas que l'itération converge.
Ensuite, l'erreur suit une récurrence d'une simplicité totale. En soustrayant de (4.21) écrite avec , et en posant :
Toute la question de la convergence tient donc dans une seule question de théorie des matrices: quand tend-il vers la matrice nulle? Ce sera le théorème 4.5. Auparavant, choisissons .
On écrit traditionnellement
où est la diagonale de , sa partie strictement triangulaire inférieure et sa partie strictement triangulaire supérieure. Les deux décompositions classiques sont alors immédiates.
On prend et . Inverser est gratuit — ce sont divisions — pourvu qu'aucun coefficient diagonal ne soit nul.
La lecture de (4.24) est parlante: on résout la -ième équation par rapport à sa -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] for j in range(n)) for i in range(n)]
if max(abs(v) for v in r) <= tol * nb:
return x, k
y = [0.0] * n
for i in range(n):
s = sum(A[i][j] * x[j] for j in range(n) if j != i)
y[i] = (b[i] - s) / A[i][i]
x = y
return x, kmaxUne itération coûte multiplications et autant d'additions pour une matrice pleine, soit environ opérations, plus pour le résidu si on le recalcule à chaque tour. Sur une matrice creuse à coefficients non nuls par ligne, ce coût tombe à . À comparer aux de la factorisation : on peut se permettre de l'ordre de itérations avant d'avoir dépensé autant qu'une élimination.
Le calcul de par (4.24) utilise 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 avec triangulaire inférieure coûte 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] for j in range(n)) for i in range(n)]
if max(abs(v) for v in r) <= tol * nb:
return x, k
for i in range(n):
s = sum(A[i][j] * x[j] for j in range(n) if j != i)
x[i] = (b[i] - s) / A[i][i] # ecriture en place: voila tout
return x, kmaxL'écriture en place x[i] = ... fait que la somme de droite lit les valeurs déjà mises à jour pour et les anciennes pour : 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.
Implémentez jacobi(A, b, tol, kmax) qui renvoie le couple formé de la solution approchée et du nombre d'itérations effectuées, avec un critère d'arrêt sur le résidu relatif en norme infinie. Le programme doit afficher 17, puis la solution du système de l'exemple 4.5 à huit décimales.
Écrivez maintenant gauss_seidel avec le même critère d'arrêt, et comparez. Le programme affiche le nombre d'itérations de Jacobi, puis celui de Gauss–Seidel, puis le rapport des deux à deux décimales, pour une tolérance de .
La récurrence (4.22) réduit la question à l'étude des puissances de .
Démonstration. Nous démontrons et , qui suffisent aux applications, et nous admettons .
: si pour une norme induite, la sous-multiplicativité (4.5) donne , donc par (4.22), et ceci quel que soit , c'est-à-dire quel que soit .
: supposons et soit une valeur propre de module maximal, un vecteur propre associé (complexe s'il le faut). Alors , 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 .
est le théorème de Householder: pour tout il existe une norme induite telle que . Sa démonstration passe par la forme normale de Jordan de 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 (formule de Gelfand, également admise): gagner un facteur demande tel que , soit .
Le critère du théorème 4.5 est exact mais coûteux: calculer 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 , 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 a pour coefficients
Sa -ième somme de ligne vaut donc, en modules,
La dominance stricte (4.26) dit exactement que pour tout . Par la formule (4.6) du théorème 4.2,
le maximum d'un nombre fini de quantités strictement inférieures à étant lui-même strictement inférieur à — c'est ici, et seulement ici, que la finitude de la dimension intervient.
L'implication du théorème 4.5 s'applique alors avec la norme : la récurrence de l'erreur (4.22) et la sous-multiplicativité donnent
ce qui est (4.27) et prouve la convergence pour tout .
Il reste l'inversibilité de . Si avec , alors , donc , d'où avec : impossible pour . Le noyau de est donc réduit à zéro et est inversible.
Une itération vérifie et . Que peut-on en conclure?
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 , avec le même que pour Jacobi (exercice 4.5). 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 appartient à la réunion des disques de centre et de rayon . 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 dominance diagonale stricte est la condition qui enferme les disques dans le disque unité.
Cet exercice est fait pour échouer d'abord. Le programme fourni lance l'itération de Jacobi sur le système témoin du cours, dont la matrice d'itération a un rayon spectral de : la méthode diverge, le critère d'arrêt n'est jamais satisfait — les composantes deviennent infinies puis indéterminées, et toute comparaison devient fausse — et l'exécution sera interrompue au bout de dix secondes. Lancez-le pour le voir, puis ajoutez le garde-fou: jacobi_surveille doit renvoyer un triplet formé de la solution, du nombre d'itérations et d'un booléen de convergence, et ne jamais boucler indéfiniment.
Gauss–Seidel calcule, pour chaque composante, une valeur censée améliorer . 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: . Avec on dépasse 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 est triangulaire inférieure de diagonale , et est triangulaire supérieure de diagonale . Donc
Le déterminant étant le produit des valeurs propres, , et le plus grand module est au moins la moyenne géométrique: . 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 un peu longue, le second sur l'analyse de la relation entre les valeurs propres de et celles de 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 avec tridiagonale d'ordre , de diagonale et de sous-diagonales — 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 à , confirme tout:
| 1,0 | 1,2 | 1,4 | 1,5 | 1,5604 | 1,6 | 1,7 | 1,8 | 1,9 | |
|---|---|---|---|---|---|---|---|---|---|
| itérations | 197 | 129 | 79 | 56 | 38 | 42 | 54 | 87 | 181 |
| 0,921 | 0,880 | 0,806 | 0,728 | 0,560 | 0,600 | 0,700 | 0,800 | 0,900 |
Trois observations, dont la troisième est celle qui sert.
D'abord, le gain est considérable: itérations contre , pour un coût par itération identique et une ligne de code de différence. Contre Jacobi, qui demande balayages sur le même système, le facteur est de .
Ensuite, la courbe est très pointue au voisinage de l'optimum: de à , soit un écart de , le nombre d'itérations passe de à . Choisir 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, coûte itérations et en coûterait environ . Pour on a exactement , une droite de pente ; en dessous, décroît beaucoup plus vite. En cas de doute, surestimez , jamais l'inverse.
Une itération a pour rayon spectral . Combien d'itérations faut-il, asymptotiquement, pour diviser l'erreur par ? Donnez le nombre d'itérations, arrondi à l'unité supérieure.
Comparons honnêtement, en opérations, sur une matrice pleine d'ordre :
L'itératif l'emporte si , c'est-à-dire si . Pour , 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. Sur une matrice pleine, l'intérêt des méthodes itératives est marginal.
Tout change sur une matrice creuse. Soit le nombre moyen de coefficients non nuls par ligne ( pour une tridiagonale, pour le laplacien 2D à cinq points, en 3D). Alors:
Chiffrons, avec une tolérance de et optimal (pour lequel le nombre d'itérations croît comme ):
| grille | Cholesky en bande | SOR optimal | mémoire: bande / creux | |
|---|---|---|---|---|
| 2 500 | op. | 149 iter., op. | / | |
| 10 000 | op. | 296 iter., op. | / | |
| 90 000 | op. | 882 iter., op. | / | |
| 1 000 000 | op. | 2 935 iter., op. | / |
Le croisement se produit vers , et l'écart s'ouvre ensuite: le coût direct croît comme , le coût itératif comme . À un million d'inconnues, SOR est fois moins cher en opérations et fois moins cher en mémoire — et c'est la mémoire qui décide en premier, car un milliard de cases, c'est 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 en secondes selon l'accès mémoire, et ce cours ne prétend pas prédire des secondes.
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 dépend des unités: multiplier par — 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 , et pas autrement. Le théorème 4.4 donne, quand (4.31) est satisfait,
Choisir sur un système de conditionnement 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 , qui ne coûte rien. Mais l'incrément et l'erreur sont reliés par
Si est proche de — 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 est énorme, et un incrément minuscule accompagne une erreur considérable. Pour , l'erreur peut valoir mille fois l'incrément. Un critère sur 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 ch4-p4.
Calculez pour la matrice dont les lignes sont, dans l’ordre, , et .
Vous pouvez afficher le corrigé directement sous chaque énoncé après avoir cherché la solution.
Soit
calculer sous forme de fraction, puis en décimal. 4. Une perturbation relative de de peut-elle rendre la solution inutilisable?
1. , , . Vérifications:
On considère le système où les lignes de sont les normales unitaires de deux droites du plan faisant entre elles l'angle :
1. , non nul pour . La formule d'inversion d'une matrice donne
On reprend avec , et .
1. , donc
On considère
1. . La dominance: , , , stricte partout. Les ratios valent , , , donc
Soit à diagonale strictement dominante en lignes, et . On note la matrice d'itération de Gauss–Seidel. On veut montrer que
1. Par définition, , ce qui s'écrit ligne par ligne
| 13,5 |
| erreur relative |
|---|
| rapport |
|---|
| 0 | 0,000000 | 0,000000 | 0,000000 | — | |
| 1 | 0,600000 | 2,000000 | −1,500000 | 0,250 | |
| 2 | 1,100000 | 1,918182 | −0,900000 | 0,200 | |
| 3 | 0,971818 | 2,018182 | −1,055682 | 0,557 | |
| 4 | 1,012955 | 1,992376 | −0,986136 | 0,249 | |
| 6 | 1,001463 | 1,999124 | −0,998202 | 0,295 | |
| 10 | 1,000020 | 1,999988 | −0,999973 | 0,334 | |
| 17 | 0,999999 | 2,000000 | −1,000000 | 0,346 |
| erreur relative |
|---|
| rapport |
|---|
| 0 | 0,000000 | 0,000000 | 0,000000 | — | |
| 1 | 0,600000 | 2,054545 | −0,879545 | 0,200 | |
| 2 | 0,981364 | 2,009256 | −0,991870 | 0,047 | |
| 3 | 0,999300 | 2,000675 | −0,999572 | 0,038 | |
| 4 | 0,999982 | 2,000037 | −0,999981 | 0,053 | |
| 5 | 1,000000 | 2,000002 | −0,999999 | 0,045 | |
| 7 | 1,000000 | 2,000000 | −1,000000 | 0,061 |
On considère le système dont la matrice a pour première ligne et pour seconde ligne , avec , et , de solution exacte . Un collègue vous présente la valeur approchée en affirmant qu'elle est excellente, «parce que le résidu est de l'ordre de ». On travaille en norme infinie et toutes les quantités sont rationnelles, donc calculables exactement.
Le déterminant vaut , donc . Les sommes de lignes de valent et ; celles de valent et .
Calculez et donnez-le en notation scientifique à trois chiffres significatifs.
Les six inégalités sont satisfaites, et la deuxième chaîne est presque saturée à droite: contre , parce que les trois composantes de ont des modules voisins, ce qui est précisément le cas d'égalité de Cauchy–Schwarz.
2. Sommes de lignes en modules: , , ; donc . Sommes de colonnes: , , ; donc . Elles coïncident parce que est symétrique: .
3. Sommes de lignes de , en modules: , , . Donc et
On vérifie au passage , et l'on remarque que ce conditionnement est tout à fait modeste: est à diagonale strictement dominante (, , ), ce qui est un bon présage — sans être une garantie, les deux notions étant distinctes.
4. Non. D'après (4.12), l'erreur relative sur la solution est majorée par , soit moins d'un demi-pour-cent: on conserve plus de deux chiffres significatifs. Un système de conditionnement inférieur à est un système sain, et c'est le cas normal.
2. Comme et sur , les deux sommes de lignes de valent : . Pour , les sommes valent et , donc, en utilisant ,
Le produit donne bien
3. Sur on a , donc et . Ainsi
fonction strictement décroissante de sur puisque la cotangente l'est. Le minimum est atteint en , où , et vaut . Deux droites perpendiculaires forment le système le mieux conditionné en norme infinie, et son conditionnement vaut 2, pas 1. (En norme 2, ce même système a : la matrice est alors orthogonale à un facteur près. La différence entre et est un artefact de la norme infinie, et elle illustre le théorème 4.1.)
4. Quand , et . Numériquement:
L'approximation donne et : elle est excellente dès le degré. Un dispositif de mesure où deux capteurs voient la scène sous un dixième de degré d'écart perd trois chiffres significatifs sur ses données, quel que soit le solveur employé.
d'où .
2. La démonstration du théorème 4.2 indique la recette: la ligne de qui réalise le maximum est la première, de coefficients . On prend donc formé des signes de cette ligne, mis à l'échelle :
Alors et : le maximum est bien atteint.
3. On a puisque : la seconde inégalité de la démonstration est elle aussi une égalité. Donc
Pour : perturbation relative , et erreur relative , soit 100 %. Une perturbation de cinq milliardièmes détruit entièrement la solution.
4. Dans l'exemple 4.3, , donc le résidu valait : le vecteur y était proportionnel à , et non à . Or la ligne de maximisante réclame les deux composantes du résidu, avec des signes opposés. Avec , seule la deuxième colonne de agit, de somme de modules au lieu de : on perd donc le facteur constaté. La borne est atteinte pour une direction de perturbation particulière, et une perturbation quelconque en réalise typiquement une fraction — ici exactement la moitié.
L'erreur est donc au moins divisée par deux par balayage, garanti.
2. Jacobi. , , ; erreur . Puis
erreur .
Gauss–Seidel. ; ; ; erreur . Puis ; ; ; erreur .
Après deux balayages, Gauss–Seidel a une erreur de contre pour Jacobi: le rapport est déjà de , et il s'établira à en régime asymptotique.
3. Pour la matrice tridiagonale d'ordre , les valeurs propres de sont , . Pour ,
La matrice étant tridiagonale donc cohéremment ordonnée, , et
4. Gagner huit chiffres demande , soit :
Le rapport Jacobi/Gauss–Seidel vaut exactement , ce qui est le cas général des matrices cohéremment ordonnées puisque . Le gain de SOR n'est ici que de supplémentaire, parce que est déjà petit: la relaxation ne paie vraiment que sur les systèmes lents, c'est-à-dire ceux de grande taille, comme le montre la figure 4.3 où le gain est d'un facteur cinq.
2. Posons et montrons par récurrence forte sur .
Initialisation, . La première somme est vide, donc
Hérédité. Supposons pour tout . Alors
Comme , on a , d'où . La récurrence est établie. (Le point fin est que la majoration utilise : c'est exactement là que la dominance stricte agit, et c'est aussi elle qui montre que la borne obtenue est meilleure que dès que .)
3. De pour tout on tire pour tout , donc par définition (4.3) de la norme induite. Le théorème 4.5 (implication ) donne alors la convergence pour tout , avec la garantie — la même que celle du théorème 4.6 pour Jacobi, donc au moins aussi bonne.
4. Non, et c'est une confusion fréquente. Ce que l'on a établi est une inégalité entre normes, , et les normes ne sont que des majorants des rayons spectraux, par (4.10). Il existe des matrices — non diagonalement dominantes — pour lesquelles Jacobi converge et Gauss–Seidel diverge, et d'autres où c'est l'inverse; le système témoin du cours en est un demi-exemple, puisque Jacobi y diverge () tandis que Gauss–Seidel y stagne ( exactement). L'égalité , qui donne le facteur 2 en nombre d'itérations, n'est démontrée que pour les matrices cohéremment ordonnées — tridiagonales par blocs, notamment. Hors de ce cadre, la comparaison se fait au cas par cas: sur le système de l'exemple 4.5, et , soit un rapport de logarithmes de et non de .