Systèmes linéaires: méthodes directes

Objectifs du chapitre

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

• expliquer pourquoi la règle de Cramer, mathématiquement irréprochable, n'est pas un algorithme, et chiffrer son coût;
• résoudre un système triangulaire par descente ou par remontée, et justifier le coût $n^2$ de chacune;
• conduire l'élimination de Gauss sur une matrice, reconnaître qu'elle produit une factorisation $LU$, et énoncer la condition sur les mineurs principaux dominants qui en garantit l'existence et l'unicité;
• établir le coût $\frac{2}{3}{n}^3$ de la factorisation et dire pourquoi on factorise une fois pour résoudre ensuite autant de seconds membres qu'on veut;
• appliquer le pivotage partiel, écrire la factorisation sous la forme , et distinguer un problème mal conditionné d'un algorithme instable sur un exemple à petit pivot;

Le problème, et une fausse solution

Tout le chapitre porte sur une seule question: étant donné $A\in {\mathbb{R}}^{n\times n}$ inversible et $\mathbf{b}\in {\mathbb{R}}^n$, calculer le vecteur $\mathbf{x}$ tel que

C'est sans doute le calcul le plus exécuté de toute l'informatique scientifique. Une simulation de structure, un circuit électrique, un écoulement, une régression, une reconstruction d'image, un pas de méthode implicite pour une équation différentielle: tous finissent par (3.1), et souvent des milliers de fois. Il vaut donc la peine de savoir exactement ce qu'il en coûte.

L'algèbre linéaire, elle, a répondu depuis longtemps: $\mathbf{x}=A^{-1}\mathbf{b}$, et la règle de Cramer donne même chaque composante en fermé,

où $A_i$ est la matrice $A$ dont la $i$-ième colonne a été remplacée par $\mathbf{b}$. Cette formule est exacte, élégante, et elle a une valeur théorique considérable — elle montre que la solution dépend de façon rationnelle des données, ce qui servira au chapitre 4. Comme algorithme, c'est une catastrophe.

Ce que $n!$ veut dire

Évaluons (3.2) honnêtement. Il faut $n+1$ déterminants d'ordre $n$. Si l'on calcule chacun par la définition de Leibniz,

la somme porte sur les $n!$ permutations et chaque terme demande $n-1$ multiplications. Un seul déterminant coûte donc $(n-1)n!$ multiplications et $n!-1$ additions.

Prenons $n=20$, une taille dérisoire — un modèle de poutre à cinq nœuds en dépasse déjà largement. On a

donc un déterminant coûte environ $4,87\cdot 10^{19}$ opérations, et les vingt et un déterminants de Cramer en réclament

Supposons — c'est une hypothèse de calcul, pas une mesure — une machine effectuant un milliard d'opérations flottantes par seconde. Il lui faudrait $1,02\cdot 10^{12}$ secondes, soit environ 32 000 ans, pour résoudre un système de vingt équations. La méthode que nous allons construire dans ce chapitre en demande $5130$: le rapport est de $2\cdot 10^{17}$.

Les systèmes triangulaires

La stratégie de tout le chapitre tient en une phrase: on ne sait résoudre facilement qu'un système triangulaire, donc on ramène tout système à des systèmes triangulaires.

Descente et remontée

Soit $L\mathbf{y}=\mathbf{b}$ avec $L$ triangulaire inférieure inversible. La première ligne s'écrit ${\ell }_{11}{y}_1=b_1$ et donne $y_1$ immédiatement. La deuxième, , ne contient plus qu'une inconnue une fois connue. De proche en proche:

C'est la substitution avant, ou descente. Symétriquement, pour $U\mathbf{x}=\mathbf{y}$ avec $U$ triangulaire supérieure inversible, on part de la dernière ligne:

C'est la substitution arrière, ou remontée. En Python, la remontée tient en six lignes:

def remontee(U, y):
    """Resout U x = y, U triangulaire superieure inversible."""
    n = len(y)
    x = [0.0] * n
    for i in range(n - 1, -1, -1):
        s = y[i]
        for j in range(i + 1, n):
            s -= U[i][j] *

Comptons. Pour l'indice $i$, la somme de (3.5) contient $n-i$ produits, donc $n-i$ multiplications et $n-i$ soustractions, suivies d'une division. Le total est

Une substitution coûte exactement $n^2$ opérations flottantes, et la descente aussi, par le même compte. Retenez ce $n^2$: il est la raison de tout ce qui suit, car il est négligeable devant le $n^3$ de la factorisation.

L'élimination de Gauss

Trois opérations qui ne changent rien

Le système (3.1) est invariant par trois opérations élémentaires sur les lignes: échanger deux lignes, multiplier une ligne par un scalaire non nul, ajouter à une ligne un multiple d'une autre. Chacune correspond à la multiplication à gauche par une matrice inversible, donc ne change ni l'ensemble des solutions, ni le rang, ni l'inversibilité. L'élimination de Gauss n'utilise que la troisième — et, à partir de la section sur le pivotage, la première.

Sur le système témoin du cours, les trois états successifs de la matrice augmentée sont ceux de la figure 3.1.

L'élimination est une factorisation déguisée

Voici l'observation qui organise tout le reste du chapitre. Pendant l'élimination, on calcule et on jette les multiplicateurs $m_{ik}$. Ne les jetons pas. Rangeons-les à la place des zéros qu'ils viennent de créer, dans le triangle inférieur, et considérons la matrice

Multiplions-la par la matrice $U$ de l'exemple 3.1. Le produit vaut

Ce n'est pas une coïncidence. L'opération $L_i\leftarrow L_i-m_{ik}{L}_k$ est la multiplication à gauche par $I-m_{ik}{\mathbf{e}}_i{\mathbf{e}}_k^{\mathsf{T}}$, dont l'inverse est — on défait la soustraction en la rajoutant. En composant toutes les étapes, on obtient avec triangulaire inférieure à diagonale unité, donc , et le miracle de bookkeeping est que est exactement la matrice (3.8): les multiplicateurs y figurent tels quels, sans changement de signe et sans mélange. C'est une propriété particulière de l'ordre dans lequel l'élimination procède, et elle rend la factorisation gratuite.

Une fois $A=LU$ obtenue, résoudre (3.1) revient à résoudre deux systèmes triangulaires:

Existence et unicité

La factorisation n'existe pas toujours. Le contre-exemple tient en quatre coefficients:

qui est inversible — c'est une matrice de permutation — et dont le premier pivot est nul: l'élimination s'arrête à la première division. Si $A=LU$, la première colonne donne ${\ell }_{11}{u}_{11}=0$ avec ${\ell }_{11}=1$, donc , puis , ce qui est impossible. La bonne condition porte sur les sous-matrices du coin supérieur gauche.

Démonstration. Par récurrence sur $n$, en construisant la factorisation bloc par bloc; l'existence et l'unicité s'obtiennent du même mouvement.

Initialisation. Pour $n=1$, $A=(a_{11})$ avec $a_{11}={\Delta }_1\ne 0$. La seule écriture avec à diagonale unité est , , et (3.10) se lit .

Hérédité. Supposons le résultat acquis pour l'ordre $k$, et écrivons $A_{k+1}$ par blocs:

Par hypothèse de récurrence, $A_k=L_k{U}_k$ de façon unique, et $\det U_k={\Delta }_k\ne 0$, donc est inversible; l'est aussi, sa diagonale valant . Cherchons et sous la forme

Cette forme est imposée: toute matrice triangulaire inférieure à diagonale unité d'ordre $k+1$ a pour bloc supérieur gauche une matrice triangulaire inférieure à diagonale unité d'ordre $k$, et de même pour $U$. En multipliant,

L'identification avec $A_{k+1}$ donne trois équations, que l'on résout dans l'ordre:

Chacune a une solution et une seule, puisque $L_k$ et $U_k$ sont inversibles. Le bloc $A_k=L_k{U}_k$ est quant à lui déjà déterminé de façon unique par l'hypothèse de récurrence. La factorisation d'ordre existe donc et est unique.

Il reste (3.10). En prenant les déterminants,

d'où $\beta =u_{k+1,k+1}={\Delta }_{k+1}/{\Delta }_k$, qui est non nul par hypothèse. Le cas $k+1=n$ donne .

L'algorithme de Doolittle

Le théorème 3.1 fournit aussi une manière de calculer $L$ et $U$ directement, en identifiant les coefficients de $A=LU$ un par un plutôt qu'en simulant l'élimination. Comme ${\ell }_{ik}=0$ pour $k>i$ et $u_{kj}=0$ pour ,

et en isolant le dernier terme on obtient les formules de Doolittle: pour $i\le j$ (partie de $U$), ${\ell }_{ii}=1$ donne

et pour $i>j$ (partie de $L$),

Il suffit de balayer dans un ordre où tout ce dont on a besoin est déjà calculé: ligne de $U$, puis colonne de $L$, et on recommence. En Python, avec la matrice écrasée sur place — $L$ sous la diagonale, $U$ sur et au-dessus, exactement comme le fait LAPACK:

def factorise_lu(A):
    """Factorisation LU sans pivotage. Renvoie (L, U) avec L a diagonale unite."""
    n = len(A)
    U = [ligne[:] for ligne in A]
    L = [[1.0 if i == j else 0.0 for j in range(n)] for i in range(n)]
    for k in range(n - 1):
        for

Ces deux boucles imbriquées sur $i$ et $j$, à l'intérieur de la boucle sur $k$, sont la raison du $n^3$ que nous allons maintenant compter exactement.

Le coût de la factorisation

Démonstration. Comptons le travail de l'étape $k$, pour $k=1,\dots ,n-1$. À cette étape, le bloc encore actif est la sous-matrice d'indices $i,j\ge k+1$, de taille $(n-k)\times (n-k)$. Pour chacune des lignes :

• une division pour former le multiplicateur $m_{ik}=a_{ik}^{(k)}/a_{kk}^{(k)}$, soit opération;

(Le coefficient d'indice $(i,k)$ lui-même n'est pas calculé: on sait qu'il vaut zéro, et on y range $m_{ik}$.) L'étape $k$ coûte donc

En posant $m=n-k$, qui décrit $n-1,n-2,\dots ,1$ quand $k$ décrit $1,\dots ,n-1$, le total vaut

Les deux sommes sont classiques:

D'où

et en réduisant au même dénominateur,

Les deux termes correctifs sont d'ordre $n^2$ et $n$, donc négligeables: $C(n)/\left(\frac{2}{3}{n}^3\right)\to 1$. $\square$

La formule (3.14) a été vérifiée en exécutant le compte brut de l'algorithme pour $n$ allant de $2$ à $1000$: les deux valeurs coïncident à l'entier près. Quelques repères, avec le rapport à l'équivalent $\frac{2}{3}{n}^3$:

Notez la dernière colonne: pour $n=1000$, une substitution coûte $10^6$ opérations contre $6,66\cdot 10^8$ pour la factorisation, soit 666 fois moins. C'est de ce rapport que vient toute la stratégie.

Une factorisation, beaucoup de seconds membres

Voici la raison d'être de $LU$, et elle n'est pas le coût brut: l'élimination de Gauss et la factorisation coûtent exactement la même chose, à quelques $n^2$ près. La différence est ailleurs.

L'élimination de Gauss agit sur la matrice augmentée: elle transforme $A$ et $\mathbf{b}$ ensemble. Si un deuxième second membre arrive plus tard — et il arrive presque toujours plus tard —, il faut tout recommencer, à $\frac{2}{3}{n}^3$ chaque fois. La factorisation, elle, ne touche pas à $\mathbf{b}$: elle produit $L$ et $U$, qui ne dépendent que de $A$, et chaque nouveau second membre ne coûte plus que les deux substitutions de (3.9), soit opérations.

Le gain se chiffre. Pour $n=1000$ et $m$ seconds membres:

Pour un seul second membre, les deux approches se valent — l'écart de $0,3\%$ est le prix des substitutions. À partir de dix, la factorisation gagne presque un facteur $m$, et à cent elle est soixante-dix-sept fois plus rapide. Or les situations à seconds membres multiples sont la règle: un pas de temps implicite résout le même système à chaque instant, une méthode de Newton pour un système non linéaire réutilise la même jacobienne plusieurs itérations de suite, une analyse de sensibilité rejoue le même modèle sous vingt chargements.

Deux sous-produits gratuits

Le déterminant. D'après (3.10), $\det A=u_{11}{u}_{22}\cdots u_{nn}$ — au signe près s'il y a eu des échanges de lignes, ce que la section suivante précisera. C'est le seul moyen raisonnable de calculer un déterminant d'ordre élevé, et il coûte $\frac{2}{3}{n}^3$ et non .

L'inverse. La $j$-ième colonne de $A^{-1}$ est la solution de $A\mathbf{z}={\mathbf{e}}_j$. On l'obtient donc par une factorisation et $n$ couples de substitutions, soit $\frac{2}{3}{n}^3+n\cdot 2n^2=\frac{8}{3}{n}^3$ opérations: Ce qui donne la règle d'hygiène la plus importante de ce chapitre.

Le pivotage partiel

Un exemple minuscule qui détruit tout

Le théorème 3.1 dit que l'élimination échoue quand un pivot s'annule. En virgule flottante, elle échoue bien avant: il suffit qu'un pivot soit petit devant les coefficients qu'il va servir à diviser. Voici le système, en deux inconnues:

Sa solution exacte est

c'est-à-dire $(1,1)$ à $10^{-17}$ près. Et le problème est excellemment conditionné: la matrice vaut $A=\left(\begin{array}{cc}\epsilon & 1\\ 1& 1\end{array}\right)$, son inverse , d'où

Un conditionnement de 4 signifie qu'une perturbation relative de $10^{-16}$ sur les données ne peut produire qu'une perturbation relative de $4\cdot 10^{-16}$ sur la solution: il n'y a aucune excuse à perdre le moindre chiffre.

Éliminons quand même sans échanger les lignes. Le multiplicateur vaut $m_{21}=1/\epsilon =10^{17}$, puis

En double précision, $10^{17}$ dépasse $2^{53}\approx 9,01\cdot 10^{15}$: l'espacement local des nombres machine y vaut 16, et $1-10^{17}$ comme s'arrondissent tous deux au nombre, . Le et le ont été purement et simplement absorbés — c'est l'annulation du chapitre 1, dans l'autre sens: une petite quantité disparaît devant une grande. La remontée donne alors

On obtient $\hat{\mathbf{x}}=(0,1)$ au lieu de $(1,1)$: une erreur relative de 100 % sur la première composante. Pas 10 %, pas trois chiffres perdus: la totalité.

Échangeons maintenant les deux lignes avant d'éliminer:

Le multiplicateur vaut cette fois $m_{21}=\epsilon =10^{-17}$, et

deux quantités qui s'arrondissent à $1$ — mais l'erreur commise est ici de $10^{-17}$ en valeur absolue sur un nombre d'ordre $1$, c'est-à-dire au niveau du bruit d'arrondi, et non $1$ sur un nombre d'ordre $10^{-17}$. On obtient $x_2=1$ puis : la solution correcte, à la précision machine.

La stratégie du pivot partiel

La recherche du maximum coûte $n-k$ comparaisons à l'étape $k$, soit $n(n-1)/2$ comparaisons en tout: un terme en $n^2$, donc gratuit devant le $n^3$ de la factorisation. C'est le meilleur rapport qualité-prix de tout le chapitre: pour un coût d'ordre inférieur, on borne tous les multiplicateurs par 1 et on rend l'échec impossible.

Il existe aussi un pivotage total, qui cherche le maximum dans toute la sous-matrice active et échange lignes et colonnes. Il donne une meilleure borne théorique sur le facteur de croissance, mais il coûte $O(n^3)$ comparaisons, ce qui n'est plus négligeable, et il n'est presque jamais utilisé: en pratique, le pivotage partiel suffit. On connaît des matrices pour lesquelles il donne un facteur de croissance de $2^{n-1}$ — elles existent, elles sont construites exprès, et elles n'apparaissent pas dans les applications. C'est un des rares endroits de ce cours où la pratique est nettement meilleure que la théorie, et il faut le dire ainsi plutôt que de prétendre que la borne est bonne.

Démonstration. Par récurrence sur $n$. Pour $n=1$, prendre $P=(1)$, $L=(1)$, $U=A$, qui est inversible.

Supposons le résultat vrai à l'ordre $n-1$ et soit $A$ inversible d'ordre $n$. Sa première colonne n'est pas nulle (sinon $A$ serait singulière), donc il existe $p$ avec

Soit $P_1$ la permutation qui échange les lignes $1$ et $p$. La matrice $P_{1}A$ a pour premier coefficient ${\tilde{a}}_{11}=a_{p1}\ne 0$, et les multiplicateurs

vérifient $|m_i|\le 1$ par le choix de $p$. Posons $\mathbf{m}=(0,m_2,\dots ,m_n{)}^{\mathsf{T}}$ et , matrice de l'élimination de la première colonne. Alors

où $A^{\prime }$ est d'ordre $n-1$. Comme $\det (M_1{P}_{1}A)=\pm \det A\ne 0$ et que ce déterminant vaut , la matrice est inversible. Par hypothèse de récurrence, il existe , , d'ordre avec et .

Posons $\tilde{P}=\left(\begin{array}{cc}1& 0^{\mathsf{T}}\\ 0& P^{\prime }\end{array}\right)$. Le point technique est le suivant: ne touche pas la première ligne, donc , et par conséquent

où ${\tilde{M}}_1$ est la même matrice d'élimination, avec ses multiplicateurs simplement permutés: permuter des lignes après coup revient à permuter les multiplicateurs et à permuter avant. Leurs modules restent $\le 1$. En multipliant l'égalité précédente par $\tilde{P}$ à gauche:

Le membre de droite se factorise:

Posons enfin $P=\tilde{P}{P}_1$ et notons $\tilde{\mathbf{m}}=\tilde{P}\mathbf{m}$. Comme ${\tilde{M}}_1^{-1}=I+\tilde{\mathbf{m}}{\mathbf{e}}_1^{\mathsf{T}}$, on obtient avec

La matrice $L$ est triangulaire inférieure à diagonale unité, ses coefficients hors diagonale sont les ${\tilde{m}}_i$ et ceux de $L^{\prime }$, tous de module $\le 1$; et $U$ est triangulaire supérieure inversible puisque ${\tilde{a}}_{11}\ne 0$ et inversible.

En pratique, $P$ ne se stocke pas comme une matrice: on garde un vecteur de $n$ entiers disant quelle ligne a été choisie à chaque étape, ce qui coûte $n$ mots au lieu de $n^2$. La résolution devient alors: permuter $\mathbf{b}$, descendre, remonter.

def resoudre(A, b):
    """Resout A x = b par elimination de Gauss avec pivot partiel."""
    n = len(b)
    M = [ligne[:] for ligne in A]
    c = list(b)
    for k in range(n - 1):
        p = max(range(k, n), key=lambda i: abs(M[i][k]))
        M[k], M[p] = M[p], M[k]
        c[k], c[p] 

Le surcoût de la recherche du pivot est de $n(n-1)/2$ comparaisons, à comparer aux $\frac{2}{3}{n}^3$ opérations flottantes de la factorisation: pour $n=1000$, environ $5\cdot 10^5$ contre , soit un millième.

Le lecteur règle le système lui-même

Dans le plan, une équation linéaire à deux inconnues est une droite, et résoudre un système $2\times 2$, c'est intersecter deux droites. Cette lecture géométrique donne tout de suite l'intuition dont le chapitre 4 fera une inégalité: quand les deux droites sont franchement obliques l'une par rapport à l'autre, leur intersection est un point net, et déplacer légèrement l'une des deux le déplace légèrement. Quand elles deviennent presque parallèles, le point d'intersection glisse le long des droites à la moindre perturbation — sans que la perturbation soit grande, et sans que le déterminant soit nul.

Faites tourner la deuxième droite jusqu'à $a_{21}\approx 1$ et $a_{22}\approx 1,05$, avec $\delta =0,2$: la solution saute de plusieurs unités pour une perturbation de $0,2$ sur un second membre d'ordre 3. Le readout «amplification observée» donne le rapport des deux erreurs relatives, et il reste toujours inférieur au conditionnement affiché à côté. C'est l'invariant de cette figure, et c'est exactement l'inégalité du chapitre 4. Remarquez enfin que le déterminant, lui, n'est pas un bon indicateur: on peut le rendre aussi petit qu'on veut en multipliant une ligne par , sans rien changer aux droites ni à la solution.

Cholesky: quand la matrice est symétrique définie positive

Une classe de matrices revient si souvent — matrices de rigidité en mécanique, matrices de covariance en statistique, matrices de Gram des moindres carrés du chapitre 7, discrétisations d'opérateurs elliptiques — qu'elle mérite son propre algorithme, deux fois plus rapide et sans pivotage.

Démonstration. Condition suffisante. Supposons $A=B{B}^{\mathsf{T}}$ avec $B$ triangulaire inférieure de diagonale strictement positive. Alors $A^{\mathsf{T}}=(B{B}^{\mathsf{T}}{)}^{\mathsf{T}}=B{B}^{\mathsf{T}}=A$, et pour ,

avec égalité seulement si $B^{\mathsf{T}}\mathbf{x}=0$; or $\det B=b_{11}\cdots b_{nn}>0$, donc est inversible et . D'où .

Condition nécessaire. Supposons $A$ SDP. Chaque sous-matrice principale dominante $A_k$ est elle-même SDP: pour $\mathbf{y}\in {\mathbb{R}}^k$ non nul, le vecteur $\mathbf{x}=(\mathbf{y},0{)}^{\mathsf{T}}\in {\mathbb{R}}^n$ est non nul et

Le déterminant d'une matrice SDP est le produit de ses valeurs propres, toutes strictement positives; donc ${\Delta }_k=\det A_k>0$ pour tout $k$. Le théorème 3.1 s'applique: $A=LU$ de façon unique, avec des pivots $u_{kk}={\Delta }_k/{\Delta }_{k-1}>0$.

Posons $D=\mathrm{diag}(u_{11},\dots ,u_{nn})$ et $\tilde{U}=D^{-1}U$, qui est triangulaire supérieure . Alors . En transposant et en utilisant :

Or ${\tilde{U}}^{\mathsf{T}}$ est triangulaire inférieure à diagonale unité et $D{L}^{\mathsf{T}}$ est triangulaire supérieure: le membre de droite est donc, lui aussi, une factorisation $LU$ de $A$. L'unicité du théorème 3.1 force

Comme tous les $u_{kk}$ sont strictement positifs, $D^{1/2}=\mathrm{diag}(\sqrt{u_{11}},\dots ,\sqrt{u_{nn}})$ est bien définie, et

avec $B$ triangulaire inférieure de diagonale $\sqrt{u_{kk}}>0$. L'unicité de $B$ découle de celle de $L$ et de : si , en factorisant chaque diagonale on retrouve deux décompositions de , donc deux factorisations , qui coïncident; les racines carrées positives étant uniques, .

L'algorithme et son coût

Les formules s'obtiennent comme celles de Doolittle, en identifiant $a_{ij}={\sum }_k{b}_{ik}{b}_{jk}$ pour $i\ge j$. Le terme diagonal donne

et les termes sous-diagonaux, pour $i>j$,

from math import sqrt

def cholesky(A):
    """Facteur de Cholesky B de A symetrique definie positive: A = B B^T."""
    n = len(A)
    B = [[0.0] * n for _ in range(n)]
    for j in range(n):
        s = A[j][j] - sum(B[j][k] * B[j][k] for k in range(j))
        B[j][j] 

Le compte est le même exercice qu'au théorème 3.2. La colonne $j$ demande $2(j-1)$ opérations et une racine pour le terme diagonal, puis, pour chacune des $n-j$ lignes en dessous, $2(j-1)$ opérations et une division. En sommant sur $j$, on trouve

soit $\frac{1}{3}{n}^3$ asymptotiquement: deux fois moins que $LU$. Pour $n=1000$, le compte exact donne $333832500$ opérations contre $666166500$, soit un rapport de . On gagne aussi la mémoire: tient dans le triangle inférieur de , et la symétrie évite d'en stocker la moitié.

Pourquoi Cholesky n'a pas besoin de pivoter

C'est la propriété la plus agréable de cette factorisation, et elle se démontre en une ligne. De $a_{ii}={\sum }_{k=1}^i{b}_{ik}^2$ on tire, pour tout $k\le i$,

Aucun coefficient du facteur ne peut dépasser la racine du plus grand terme diagonal de $A$: le facteur de croissance est borné a priori, indépendamment de la matrice, donc il n'y a rien à surveiller pendant l'élimination et aucun échange de lignes à faire. De plus, les pivots $b_{jj}^2={\Delta }_j/{\Delta }_{j-1}$ sont strictement positifs et la racine carrée de (3.22) est toujours définie — en arithmétique exacte.

Matrices bandes et matrices creuses

La bande, et ce qu'elle fait gagner

Ces matrices sont omniprésentes: discrétiser une équation différentielle par différences finies donne une matrice tridiagonale (chapitre 9), un problème bidimensionnel sur une grille $m\times m$ donne une matrice de largeur $m$, une spline cubique du chapitre 7 conduit à un système tridiagonal.

La propriété décisive est que l'élimination de Gauss sans pivotage préserve la bande. Le raisonnement est immédiat sur (3.7): à l'étape $k$, seules les lignes $i\le k+p$ ont un coefficient non nul en colonne $k$, donc seules celles-là sont modifiées; et la ligne $k$ n'a de coefficients non nuls que jusqu'à la colonne $k+q$, donc la modification ne s'étend pas au-delà. Aucun zéro situé hors de la bande n'est jamais touché. Le facteur $L$ hérite de la largeur inférieure $p$ et de la largeur supérieure .

Le coût suit. À l'étape $k$, le bloc actif n'est plus $(n-k)\times (n-k)$ mais au plus $p\times q$, d'où

Pour $n=10000$ et $p=q=50$, le compte exact de l'algorithme donne $5,03\cdot 10^7$ opérations, contre $6,67\cdot 10^{11}$ pour la même matrice traitée comme pleine: . Le cas tridiagonal est encore plus favorable: la factorisation et les deux substitutions se fondent en un algorithme unique, dit de , dont le coût exact est opérations — pour , contre en plein. On est passé de à .

Le creux, et le remplissage

Une matrice creuse est une matrice dont la grande majorité des coefficients sont nuls, sans que ces zéros forment nécessairement une bande. On ne stocke alors que les coefficients non nuls, avec leurs indices, et l'on espère que la factorisation restera creuse elle aussi.

Cet espoir est souvent déçu, et c'est le phénomène central de l'algèbre linéaire creuse.

Le mécanisme est donc purement combinatoire: c'est le motif des zéros, et lui seul, qui décide du remplissage. Et il dépend violemment de l'ordre dans lequel on numérote les inconnues — qui est pourtant arbitraire, puisque renuméroter les inconnues et les équations ne change ni le système ni sa solution.

L'exemple de la figure 3.3 est extrême, et c'est pour cela qu'il est instructif. La matrice flèche

a $3n-2$ coefficients non nuls, soit 28 pour $n=10$. La première étape de l'élimination utilise la ligne 1, qui est pleine, pour modifier toutes les lignes suivantes, qui ont toutes un coefficient en colonne 1: chaque couple $(i,j)$ avec $i,j\ge 2$ reçoit une contribution, et la sous-matrice restante devient pleine d'un coup. Les facteurs ont $n^2=100$ coefficients non nuls, soit 72 de remplissage, et la factorisation coûte le prix plein.

Renumérotons maintenant les inconnues dans l'ordre inverse. La flèche pointe vers le bas à droite: chaque ligne $i<n$ n'a que deux coefficients, en $(i,i)$ et $(i,n)$. À l'étape $k$, la seule ligne à modifier est la ligne $n$, et le seul coefficient touché est en colonne $n$ — où il est déjà non nul. Le remplissage est nul, les facteurs ont exactement les 28 coefficients de départ, et le coût est linéaire en $n$.

Même matrice, mêmes équations, même solution: seul l'ordre de numérotation a changé, et le coût est passé de $O(n^3)$ à $O(n)$. C'est pourquoi tout solveur creux sérieux commence par une phase d'analyse symbolique, qui ne regarde que le motif des zéros et cherche une permutation $P$ minimisant le remplissage de $PA{P}^{\mathsf{T}}$, avant toute arithmétique. Trouver la permutation optimale est un problème NP-difficile, et l'on utilise donc des heuristiques, qui forment une petite bibliothèque classique: Cuthill–McKee inverse, qui regroupe les coefficients près de la diagonale pour fabriquer une bande étroite; le degré minimum, qui élimine à chaque étape l'inconnue créant le moins de remplissage; la dissection emboîtée, qui coupe récursivement le graphe de la matrice en deux moitiés séparées par un petit ensemble de sommets. Ces trois familles portent le nom de leur idée et non d'un gain chiffré, parce que ce gain dépend entièrement de la structure du problème: il peut être nul, ou valoir plusieurs ordres de grandeur comme dans la figure 3.3.

Ce qu'une bibliothèque appelle vraiment

Il vaut la peine de savoir ce qui se passe derrière un appel comme numpy.linalg.solve(A, b) — d'autant que ce cours vous demande d'écrire les boucles vous-même précisément pour que cet appel cesse d'être une boîte noire.

Sous presque toutes les bibliothèques d'algèbre linéaire dense se trouve LAPACK, écrite en Fortran depuis 1992, et sous LAPACK les BLAS, trois niveaux de briques élémentaires: niveau 1 pour les opérations vecteur-vecteur, niveau 2 pour matrice-vecteur, niveau 3 pour matrice-matrice. Les routines pertinentes ici portent des noms de six lettres dont la première, d, indique la double précision:

L'algorithme de dgetrf est exactement celui de ce chapitre, à une transformation près: il est par blocs. Au lieu de traiter une colonne à la fois, il factorise un panneau de $b$ colonnes (typiquement $b=64$), puis met à jour le reste de la matrice par un seul produit matriciel, qui est une opération de niveau 3. Le nombre d'opérations flottantes est inchangé, $\frac{2}{3}{n}^3$: ce qui change est le nombre d'accès à la mémoire, car un produit matriciel réutilise chaque donnée chargée $b$ fois au lieu d'une. C'est une leçon générale de calcul scientifique, et elle ne contredit pas le comptage d'opérations de ce chapitre — elle le complète: à nombre d'opérations égal, l'organisation des données décide de la vitesse réelle.

Du côté creux, les solveurs directs de référence (SuperLU, UMFPACK, MUMPS, CHOLMOD pour le cas SDP) enchaînent toujours les mêmes trois phases: analyse symbolique — la renumérotation de la section précédente, qui ne fait aucune arithmétique —, factorisation numérique, puis résolution pour chaque second membre. La séparation en trois phases est la traduction logicielle de la règle de ce chapitre: ce qui ne dépend que de $A$ se calcule une fois.

Quelques réflexes, pour finir:

• appeler solve, jamais inv — quatre fois moins cher et plus précis;
• si plusieurs seconds membres, conserver la factorisation (scipy.linalg.lu_factor puis lu_solve);
• si la matrice est SDP, le dire à la bibliothèque (cho_factor, assume_a="pos"): on gagne le facteur 2 du théorème 3.4 et la garantie de stabilité sans pivotage;
• si elle est creuse, la stocker creuse, et non comme un tableau plein de zéros: une matrice $10^5\times 10^5$ pleine demanderait 80 gigaoctets;
• et toujours, à la fin, regarder le résidu.

Synthèse

• La règle de Cramer et l'inverse explicite sont des outils de démonstration, pas des algorithmes: à $n=20$, Cramer avec la formule de Leibniz demande $1,02\cdot 10^{21}$ opérations, contre $5130$ pour la factorisation — un rapport de $2\cdot 10^{17}$. Calculer $A^{-1}$ pour résoudre un système coûte quatre fois trop cher et donne un résultat moins précis.

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 (méthodes directes pour les systèmes linéaires, factorisations $LU$ et de Cholesky, matrices bandes).
• Rappaz, J. et Picasso, M., Introduction à l'analyse numérique, Presses polytechniques et universitaires romandes, Lausanne, chap. 2.
• Trefethen, L. N. et Bau, D., Numerical Linear Algebra, SIAM, Philadelphie, leçons 20 à 23 (élimination de Gauss, pivotage, stabilité et facteur de croissance).
• Golub, G. H. et Van Loan, C. F., Matrix Computations, Johns Hopkins University Press, Baltimore, chap. 3 et 4 (algorithmes par blocs, matrices bandes, Cholesky).
• Higham, N. J., Accuracy and Stability of Numerical Algorithms, 2e éd., SIAM, Philadelphie, chap. 9 (analyse d'erreur de l'élimination de Gauss et du pivotage).
• Davis, T. A., Direct Methods for Sparse Linear Systems, SIAM, Philadelphie (renumérotation, remplissage et solveurs creux).