Résolution de systèmes linéaires

Par des méthodes directes

 

Elimination de Gauss
 

Le processus d'élimination de Gauss est probablement le plus simple que l'on puisse envisager. Le principe est d'exprimer une des inconnues en fonction des autres, de substituer la valeur obtenue dans les autres équations, de répéter l'opération jusqu'à l'obtention d'une équation à une seule inconnue, ce qui permet ensuite d'en déduire l'ensemble des autres inconnues. Essayons par exemple de résoudre le système :

x - 7y = -11
2x + 3y = 12

La première équation nous permet d'écrire x = 7y - 11, ce qui donne en réinjectant cette expression dans la deuxième équation :

2(7y-11) + 3y = 12, d'où 17y = 34 et par conséquent y = 2. On remonte ensuite facilement à x = 3.

Voici une description simplifiée de la partie "élimination" de l'algorithme :

Pour i = 1 à n-1,
    pour j = i+1 à n
        a_ij = a_ij / a_ii
    fin pour
    b_i = b_i / a_ii
    pour k = i+1 à n
        pour j = i+1 à n
            a_kj = a_kj - a_ki.a_ij
        fin pour
        b_k = b_k - a_ki.b_i
    fin pour
fin pour
b_n = b_n/a_nn

Il suffit ensuite de résoudre le système triangulaire supérieur ainsi obtenu. Le nombre d'opérations (multiplications et divisions) de cet algorithme est approximativement de (n³/3 - n/3) pour la matrice A et de n² pour chaque vecteur b.

Si toutes les sous-matrices de a ne sont pas régulières, on risque de tomber sur un pivot nul. Il faut alors intégrer une procédure de choix du plus grand pivot. Le pivotage partiel consiste à échanger deux équations dans le but d'avoir le plus grand pivot possible en valeur absolue. Le pivotage complet consiste à procéder à un échange de lignes ET de colonnes, mais il requiert un nombre d'opérations de l'ordre de n³, contre n² pour le partiel. Il n'est donc jamais pratiqué.

Lorsque deux systèmes linéaires assez voisins aboutissent à des résultats très différents, on parle de mauvais conditionnement. Le conditionnement est donc la mesure de la sensibilité d'un système à de faibles modifications de l'un de ses éléments (que ce soit la matrice ou le vecteur du second membre). Supposons que l'on parte de A.x = b. Puisque ||AB|| <,= ||A||.||B||, on a ||b|| <,= ||A||.||x||. Altérons légèrement b en b+db. La solution x+dx du nouveau système est telle que A.(x+dx) = b+db. En soustrayant Ax = b à cette dernière équation, on a dx = A^-1.db, donc  ||dx|| <,= ||A^-1||.||db||. En multipliant ceci par ||b|| <,= ||A||.||x||, on obtient ||b||.||dx|| <,= ||A||.||x||.||A^-1||.||db||. En divisant par ||b||.||x||, on obtient :

(||dx|| / ||x||) <,= ||A||.||A^-1||.(||db|| / ||b||)

qui montre que "l'erreur" relative sur la solution x est bornée par le produit de l'erreur relative sur le vecteur b et de la quantité ||A||.||A^-1||.
Une mesure (parmi d'autres) du conditionnement est donnée par r = |||A|||.|||A^-1|||, où |||A||| est ici la norme spectrale (i.e la plus petite valeur propre) de A. Le rayon spectral r(A) est la valeur absolue du rapport entre la plus grande et la plus petite valeur propre d'une matrice A.
La règle suivante peut être retenue : sur un calculateur à p chiffres significatifs, le nombre de chiffres significatifs de la solution est E[p-log₁₀ r(A)], où E[] est la partie entière.
On comprend en tout cas aisément que plus un système est mal conditionné, plus la solution sera sensible aux erreurs d'arrondi. Ceci n'implique par contre évidemment pas qu'un système bien conditionné soit systématiquement peu sensible aux erreurs d'arrondi ...

Système surdéterminé

Dans le cas d'un système surdéterminé, c'est-à-dire du type A.x = b avec A matrice [m,n] et m > n, on cherche à déterminer le vecteur x qui minimise la norme de A.x - b. En pré-multipliant A.x = b par A^t , on obtient A^t.A.x = A^t.b, où A^t.A est une matrice [n,n] symétrique définie positive. Il est facile de montrer que la solution de ce système est également solution du problème aux moindres carrés initial.

Décomposition LU

Si toutes les sous-matrices de A sont régulières, A est décomposable en un produit L.U, où L et U sont des matrices triangulaires inférieure et supérieure, l'une d'entre elles (on prendra U par la suite, selon l'algorithme proposé par Crout) contenant des 1 sur sa diagonale. U correspond à la matrice obtenue par le processus d'élimination de Gauss.

La forme particulière de ces matrices permet d'établir un algorithme où le stockage de U et L se fera directement "dans" A.

Pour i = 2 à n,
    a_1i = a_1i/a₁₁
fin pour
pour k = 2 à n-1
    a_kk = a_kk - Σ_j=1->k-1 a_kj. a_jk
    pour i = k+1 à n
        a_ik = a_ik - Σ_j=1->k-1 a_ij.a_jk
        a_ki = (a_ki - Σ_j=1->k-1 a_kj.a_ji) / a_kk
    fin pour
fin pour
a_nn = a_nn - Σ_j=1->n-1 a_nj.a_jn  

Lorsque l'on souhaite résoudre un système A.x = b, la décomposition LU permet de résoudre successivement L.y = b puis U.x = y. La résolution du système linéaire pour des seconds membres différents est donc beaucoup plus rapide.

Le nombre d'opérations pour la décomposition LU d'une matrice n x n est de (n³ - n)/3 multiplications/divisions, (2n³ - 3n² + n)/6 additions/soustractions, et les remontées/descentes coûtent n² mult./div. et (n² -n) add./soustr. On peut donc dire que le gros du travail est constitué par la décomposition, qui évolue en n³/3.

Si l'on souhaite inverser une matrice, il est souvent préférable de passer par une décomposition LU pour résoudre n systèmes linéaires (avec un coût total en 4n³/3) plutôt que de calculer réellement l'inverse.

Décomposition de Cholesky

Rappelons tout d'abord que toutes les sous-matrices d'une matrice symétrique définie positive le sont également, et sont donc régulières.

Le théorème de Cholesky permet de dire que toute matrice carrée A symétrique définie positive admet une décomposition unique A = L.L^t où L est une matrice triangulaire inférieure à valeurs diagonales positives.

L'algorithme peut s'écrire de la manière suivante :

a₁₁ = a₁₁^1/2
Pour i = 2 à n,
    a_i1 = a_i1/a₁₁
fin pour
pour k = 2 à n-1
    a_kk = (a_kk - Σ_j=1->k-1 a_kj²)^1/2
    pour i = k+1 à n
        a_ik = (a_ik - Σ_j=1->k-1 a_ij.a_kj) / a_kk
    fin pour
fin pour
a_nn = (a_nn - Σ_j=1->n-1 a_nj²)^1/2
 

Matrices-bande

Si on note l la demi-largeur de bande d'une matrice, on peut montrer que le nombre d'opérations nécessaires pour une décomposition LU (ou LL^t) évolue comme l²n. Pour les matrices tridiagonales, qui apparaissent notamment dans la résolution d'équations différentielles par des méthodes implicites, l'algorithme de Thomas qui dérive du procédé de Gauss permet de résoudre le système en 5n-4 multiplications (et seulement 3n-2 opérations par vecteur b supplémentaire). Le pivotage ne peut être utilisé avec l'algorithme de Thomas parce qu'il détruit la tridiagonalité, mais comme les systèmes tridiagonaux sont généralement à diagonale dominante le pivotage n'est pas nécessaire.

 

Par des méthodes itératives

 

Méthodes de Jacobi et Gauss-Seidel

Dans de nombreux problèmes numériques la matrice A est dite creuse, ce qui signifie qu'elle contient une majorité de termes nuls. Dans ces cas là les méthodes directes de résolution trouvent fréquemment des adversaires coriaces (en termes de rapidité) dans la famille des méthodes itératives. Celles-ci consistent à partir d'un vecteur initial (plus ou moins bien choisi) et à construire une suite de vecteurs qui converge vers la solution du système linéaire. On décompose pour ce faire A sous la forme A = D - E - F, où D est la matrice diagonale formée des a_ii, -E la matrice triangulaire inférieure contenant les éléments sous-diagonaux de A, et -F la matrice triangulaire supérieure contenant les éléments sur-diagonaux de A.

A.x = b se réécrit alors D.x = b + (E+F).x, système diagonal, ce qui permet d'obtenir la méthode de Jacobi :

x_i = (b_i - Σ_j=1->n,j#i a_ij.x_j) / a_ii

Petit exemple :

4x₁ - x₂           = 100
-x₁ + 4x₂ - x₃ = 100 
        -x₂ + 4x₃ = 100

Ce système peut se réécrire sous la forme :

4x₁ = x₂ + 100
4x₂ = x₁ + x₃ + 100
4x₃ = x₂ + 100 

En partant d'un vecteur initial (0, 0, 0), on peut donc itérer :

4x₁ = 100           x₁ = 25
4x₂ = 100   =>    x₂ = 25
4x₃ = 100           x₃ = 25

puis :

4x₁ = 125           x₁ = 31.25
4x₂ = 150   =>    x₂ = 37.5
4x₃ = 125           x₃ = 31.25

etc ...

La méthode de Gauss-Seidel consiste à résoudre (D-E).x = b + F.x, qui est un système triangulaire inférieur. Si l'on note x⁽ⁿ⁾ l'itéré de rang n, l'écriture composante par composante montre que cela revient à prendre en compte les termes de x⁽ⁿ⁺¹⁾ venant d'être calculés pour obtenir les composantes restantes. Voilà pourquoi on parle parfois de méthode des "itérations successives", par rapport au nom de "méthode des itérations simultanées" qui désigne parfois celle de Jacobi.

Les méthodes itératives définies par une relation du type x⁽ⁿ⁺¹⁾ = M.x⁽ⁿ⁾ + p sont convergentes si le rayon spectral r(M) = max_j |l_j|, où l_j est la j-ième valeur propre de A est strictement inférieur à 1. On peut montrer que :
- si A est symétrique définie positive, la méthode de Gauss-Seidel est convergente.
- si A et 2D - A sont symétriques définies positives, la méthode de Jacobi est convergente.
- si A est à diagonale dominante, les méthodes de Gauss-Seidel, Jacobi, et SOR (voir ci-dessous) convergent. En revanche, si A ne peut être transformée en une matrice à diagonale dominante, le système peut converger pour certains vecteurs initiaux mais la convergence n'est pas garantie quel que soit ce vecteur.

La méthode de Gauss-Seidel converge généralement (mais pas toujours) plus vite que celle de Jacobi.

Les erreurs d'arrondi à une itération donnée ne dépendent pas des erreurs d'arrondi à l'itération précédente, alors que pour les méthodes directes on a vu que les erreurs d'arrondi polluent progressivement la solution.

Méthodes de relaxation

En utilisant les méthodes qui viennent d'être présentées on se rend compte que la correction apportée au vecteur solution à chaque itération a tendance à être sous-estimée. En d'autres termes, le vecteur ne converge pas "assez vite" vers la solution. De là vient l'idée de "corriger la correction", à l'aide d'un facteur multiplicatif
En reprenant la décomposition A = D - E - F, et si on note w un réel non nul que l'on appellera paramètre de relaxation, on peut écrire

(D/w - E).x⁽ⁿ⁺¹⁾ = ((1-w).D/w + F).x⁽ⁿ⁾ + b

qui nécessite comme Gauss-Seidel la résolution d'un système triangulaire à chaque itération. On parle de sous- ou de surrelaxation selon que w<1 ou w>1. On peut montrer que cette méthode converge pour 0 < w < 2, et que le paramètre optimal est donné par

w_opt = 2 / (1 + sqrt(1-r²(J)))

où r(J) est le rayon spectral de la matrice J = D^-1(E+F) de la méthode de Jacobi.

La surrelaxation successive symétrique (SSOR) consiste à faire jouer le même rôle aux matrices E et F, en introduisant un vecteur intermédiaire y⁽ⁿ⁾ :

(D/w - E).y⁽ⁿ⁾ = ((1-w).D/w + F).x⁽ⁿ⁾ + b
(D/w - F).x⁽ⁿ⁺¹⁾ = ((1-w).D/w + E).y⁽ⁿ⁾ + b

Méthodes du gradient et du gradient conjugué

On définit l'application L(x) = ½.x^t.A.x - b^t.x. Si A est symétrique définie positive et si x est solution de A.x = b, alors x minimise L.
Si x_n est le n-ième itéré, on peut chercher l'itéré suivant sous la forme

x_n+1 = x_n + q_n+1.w_n+1

où q_n+1 est un réel est w_n+1 un vecteur donné. On cherche le réel q_n+1 qui permet de minimiser L(x_n+1). On peut montrer qu'il est donné par

q_n+1 = (w_n+1)^t.(b - A.x_n) / [(w_n+1)^t.A.w_n+1]

Il faut en outre savoir quelle direction de descente w choisir. On peut choisir la direction de plus grande pente, qui correspond au gradient de L. On constate de suite que grad(L(x_n)) = A.x_n - b. Le vecteur w est donc égal au résidu. Si A est symétrique définie positive, cette méthode converge en théorie. Mais la convergence peut être très lente, et elle peut être perdue à cause des erreurs d'arrondis.

La solution consiste à corriger la direction de descente en intégrant une composante fonction de la direction à l'itération précédente :

w_n+1 = -r_n + b_nw_n

où b_n est un réel calculé de manière à minimiser ||x - x_n||. On montre que

b_n = (r_n)^t.A.w_n / [(w_n)^t.A.w_n]

Le résidu ainsi obtenu à l'itération n est orthogonal aux n résidus précédents, ce qui permet de prouver que si A est une matrice NxN symétrique définie positive, la méthode du gradient conjugué converge en au plus N itérations.

Cette méthode est cependant fréquemment mal conditionnée, la vitesse de convergence dépendant directement du rayon spectral de A. On prémultiplie donc le système par l'inverse d'une matrice C, de façon à ce que le rayon spectral de C^-1A soit plus petit que celui de A. C^-1A n'étant pas symétrique, il est nécessaire de procéder à une décomposition de Cholesky de C, ce qui permet d'écrire

L^-tL^-1A.x = L^-tL^-1.b

puis

L^-1AL^-tL^t.x = L^-1.b

ce qui revient à un système A*.x* = b*, en notant L^-1AL^-t = A^*, L^t.x = x^*, et L^-1.b = b^*, et où L^-1AL^-t est symétrique définie positive. C'est le choix de la matrice de préconditionnement C qui va ... conditionner la rapidité de convergence de l'algorithme.

Il existe également des méthodes proches du gradient conjugué et qui soient adaptées aux matrices non symétriques : cf. GMRES et BICG-Stab
 
Critères de terminaison d'une méthode itérative

L'éventuelle convergence d'une méthode se juge à la variation entre deux itérés successifs. On peut formuler un critère en termes de variation relative ou de variation absolue. A moins de connaître l'ordre de grandeur approximatif de la solution, le premier choix est bien évidemment préférable.