Changement de pivot dans la méthode de Gauss

Nous avons déjà signalé que la méthode de Gauss échoue si un pivot s’annule. Voir l’exemple.

Dans ce cas, on doit recourir à une technique dite de changement de pivot qui consiste à échanger des lignes (ou des colonnes) du système de manière à ce qu’aucun pivot ne soit nul.

Exemple: Échec de la méthode de Gauss

Revenons à la matrice chap3/gauss.adoc#eq19 pour laquelle la méthode de Gauss donne un pivot nul à la seconde étape. En échangeant simplement la deuxième et la troisième ligne, on trouve un pivot non nul et on peut exécuter une étape de plus. Le système obtenu est équivalent au système de départ, et on constate qu’il est déjà triangulaire supérieur. En effet,

\[\mathrm{A}^{(2)}=\left[\begin{array}{rrr} 1 & 2 & 3 \\ 0 & -6 & -12 \\ 0 & 0 & -1 \end{array}\right]=\mathrm{U}\]

et les matrices de transformation sont données par

\[\mathrm{M}_1=\left[\begin{array}{ccc} 1 & 0 & 0 \\ -2 & 1 & 0 \\ -7 & 0 & 1 \end{array}\right], \quad \mathrm{M}_2=\left[\begin{array}{lll} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array}\right]\]

D’un point de vue algébrique, ayant effectué une permutation des lignes de \(A\), l’égalité \(A=M_1^{-1} M_2^{-1} U\) doit être remplacée par \(A=M_1^{-1} P M_2^{-1} U\) où \(P\) est la matrice de permutation

\[\mathrm{P}=\left[\begin{array}{lll} 1 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0 \end{array}\right]\]
Changement de pivot partiel

La stratégie de pivot adoptée ci-dessus peut être généralisée en recherchant, à chaque étape \(k\) de l’élimination, un pivot non nul parmi les termes de la sous-colonne \(\mathrm{A}^{(k)}(k: n, k)\). On dit alors qu’on effectue un changement de pivot partiel (par ligne).

On peut voir à partir de xref:chap3/gauss.adoc#eq17[la définition des \(m_{i k}\) qu’une grande valeur de \(m_{i k}\) (provenant par exemple d’un petit pivot \(\left.a_{k k}^{(k)}\right)\) peut amplifier les erreurs d’arrondi affectant les termes \(a_{k j}^{(k)}\).

Par conséquent, afin d’assurer une meilleure stabilité, on choisit comme pivot l’élément de la colonne \(\mathrm{A}^{(k)}(k: n, k)\) le plus grand en module et le changement de pivot partiel est généralement effectué à chaque étape, même si ce n’est pas strictement nécessaire (c’est-à-dire même s’il n’y a pas de pivot nul).

Changement de pivot total

Une méthode alternative consiste à rechercher le pivot dans l’ensemble de la sous-matrice \(\mathrm{A}^{(k)}(k\) : \(n, k: n\) ), effectuant alors un changement de pivot total

le changement de pivot partiel ne requiert qu’un surcoût d’environ \(n^2\) tests, alors que le changement de pivot total en nécessite environ \(2 n^3 / 3\), ce qui augmente considérablement le coût de la méthode de Gauss.

Analysons comment la stratégie de pivot partiel affecte la factorisation LU induite par la méthode de Gauss. Lors de la première étape de l’élimination de Gauss avec changement de pivot partiel, après avoir trouvé l’élément \(a_{r 1}\) de module maximum dans la première colonne, on construit la matrice de permutation élémentaire \(\mathrm{P}_1\) qui échange la première et la \(r\)-ième ligne (si \(r=1, \mathrm{P}_1\) est la matrice identité). On crée ensuite la première matrice de transformation de Gauss \(\mathrm{M}_1\) et on pose \(A^{(2)}=M_1 P_1 A^{(1)}\). On procède de même avec \(A^{(2)}\), en cherchant les nouvelles matrices \(P_2\) et \(\mathrm{M}_2\) telles que

\[\mathrm{A}^{(3)}=\mathrm{M}_2 \mathrm{P}_2 \mathrm{~A}^{(2)}=\mathrm{M}_2 \mathrm{P}_2 \mathrm{M}_1 \mathrm{P}_1 \mathrm{~A}^{(1)}\]

Après avoir effectué toutes les étapes de l’élimination, la matrice triangulaire supérieure \(U\) obtenue est donnée par

\[\mathrm{U}=\mathrm{A}^{(n)}=\mathrm{M}_{n-1} \mathrm{P}_{n-1} \ldots \mathrm{M}_1 \mathrm{P}_1 \mathrm{~A}^{(1)} .\]

En posant \(\mathrm{M}=\mathrm{M}_{n-1} \mathrm{P}_{n-1} \ldots \mathrm{M}_1 \mathrm{P}_1\) et \(\mathrm{P}=\mathrm{P}_{n-1} \ldots \mathrm{P}_1\), on a \(\mathrm{U}=\mathrm{MA}\) et donc \(\mathrm{U}=\left(\mathrm{MP}^{-1}\right) \mathrm{PA}\). On vérifie facilement que la matrice \(\mathrm{L}=\mathrm{PM}^{-1}\) est triangulaire inférieure (avec des 1 sur la diagonale), de sorte que la factorisation LU s’écrit

\[\mathrm{PA}=\mathrm{LU}\]
On ne doit pas s’inquiéter de la présence de l’inverse de \(\mathrm{M}\), car \(\mathrm{M}^{-1}=\mathrm{P}_1^{-1} \mathrm{M}_1^{-1} \ldots \mathrm{P}_{n-1}^{-1} \mathrm{M}_{n-1}^{-1}\), \(\mathrm{P}_i^{-1}=\mathrm{P}_i^T\) et \(\mathrm{M}_i^{-1}=2 \mathrm{I}_n-\mathrm{M}_i\).

Une fois calculées les matrices \(\mathrm{L}, \mathrm{U}\) et \(\mathrm{P}\), la résolution du système initial se ramène à la résolution des systèmes triangulaires \(\mathrm{Ly}=\mathrm{Pb}\) et \(\mathrm{Ux}=\mathbf{y}\).

Remarquer que les coefficients de la matrice \(L\) coïncident avec les multiplicateurs calculés par une factorisation \(LU\) de la matrice \(PA\) sans changement de pivot.

Si on adopte une stratégie de pivot total, à la première étape, une fois trouvé l’élément \(a_{q r}\) de plus grand module dans la sous-matrice \(\mathrm{A}(1: n, 1: n)\), on doit échanger la première ligne et la première colonne avec la \(q\)-ième ligne et la \(r\)-ième colonne. Ceci conduit à la matrice \(\mathrm{P}_1 \mathrm{~A}^{(1)} \mathrm{Q}_1\), où \(\mathrm{P}_1\) et \(\mathrm{Q}_1\) sont respectivement des matrices de permutation de lignes et de colonnes. A ce stade, la matrice \(M_1\) est donc telle que \(A^{(2)}=M_1 P_1 A^{(1)} Q_1\). En répétant ce processus, on obtient à la dernière étape

\[\mathrm{U}=\mathrm{A}^{(n)}=\mathrm{M}_{n-1} \mathrm{P}_{n-1} \ldots \mathrm{M}_1 \mathrm{P}_1 \mathrm{~A}^{(1)} \mathrm{Q}_1 \ldots \mathrm{Q}_{n-1},\]

au lieu de [eq24].

Dans le cas d’un changement de pivot total, la factorisation \(LU\) devient

\[\mathrm{PAQ}=\mathrm{LU}\]

où \(\mathrm{Q}=\mathrm{Q}_1 \ldots \mathrm{Q}_{n-1}\) est une matrice de permutation prenant en compte toutes les permutations effectuées. Par construction, la matrice \(L\) est encore triangulaire inférieure, et ses éléments ont tous un module inférieur ou égal à 1 . Comme pour le changement de pivot partiel, les éléments de L sont les multiplicateurs produits par la factorisation LU de la matrice PAQ sans changement de pivot. Une fois calculées les matrices \(\mathrm{L}, \mathrm{U}, \mathrm{P}\) et \(\mathrm{Q}\), pour résoudre le système linéaire on note qu' on peut écrire

\[\mathrm{Ax}=\mathrm{b} \Leftrightarrow \underbrace{\mathrm{PAQ}}_{\mathrm{LU}} \underbrace{\mathrm{Q}^{-1} \mathbf{x}}_{\mathbf{x}^*}=\mathrm{Pb} .\]

Alors, on se ramème à la résolution des systèmes triangulaires \(\mathrm{Ly}=\mathrm{Pb}, \mathrm{Ux}^*=\mathrm{y}\) et finalement on trouve \(\mathbf{x}=\mathbf{Q} \mathbf{x}^*\).

1. Python

Dans Python, en utilisant la bibliothèque numpy, on peut calculer la factorisation LU d’une matrice A comme suit :

import numpy as np
P, L, U = np.linalg.lu(A)

La matrice P est la matrice de permutation. Lorsque la matrice P coïncide avec la matrice identité, les matrices L et U sont les matrices recherchées (telles que \(LU = A\)). Autrement, on a \(LU = PA\).

Exemple: Factorisation LU en Python

Considérons la matrice suivante :

\[A = \begin{bmatrix} 4 & 3 & 2 & 1 \\ 6 & 3 & 4 & 5 \\ 8 & 7 & 6 & 5 \\ 2 & 1 & 2 & 3 \end{bmatrix}\]

Appliquons la factorisation LU en Python :

import numpy as np
from scipy.linalg import lu

A = np.array([[4, 3, 2, 1], [6, 3, 4, 5], [8, 7, 6, 5], [2, 1, 2, 3]])
P, L, U = lu(A)

print(f"P: {P}")
print(f"L: {L}")
print(f"U: {U}")

# Vérification que LU = PA
assert np.allclose(np.dot(L, U), np.dot(P, A)), "Erreur dans la factorisation!"
print("Vérification réussie: LU = PA!")

Avec cette exécution, nous aurons les matrices L, U et P imprimées à l’écran. Si la factorisation est correcte, le message "Vérification réussie: LU = PA!" sera également affiché.

Mesurons la performance (en temps) de la méthode LU avec pivot partiel en augmentant la taille de la matrix et en utilisant les matrices du Laplacian

import numpy as np
from scipy.linalg import lu
import time

# Tailles des matrices
Ns = [10, 100,  1000, 2000] # 3000, 4000, 5000, 6000, 7000]
Ts = []
nsamples = 4
for N in Ns:
    h = float(1./N)
    x = np.linspace(start=h, stop=1 - h, num=N - 1)
    # Construct the matrix A
    A = (N**2) * (np.diag(2 * np.ones(N - 1)) - np.diag(np.ones(N - 2), 1) - np.diag(np.ones(N - 2), -1))

    tavg = 0
    for i in range(nsamples):
        # Mesure du temps de la factorisation LU
        start = time.time()
        P, L, U = lu(A)
        end = time.time()
        tavg += end - start

    tavg/=nsamples
    # Affichage du temps de calcul
    Ts.append( tavg )
    print(f"N = {N:4d} | Temps de calcul: {tavg:.4f} s")
Results 
N =   10 | Temps de calcul: 0.0001 s
N =  100 | Temps de calcul: 0.0003 s
N = 1000 | Temps de calcul: 0.0280 s
N = 2000 | Temps de calcul: 0.1781 s

Nous tracons le graphe en loglog des résultats avec plotly

Tracé LogLog du temps de calcul de la factorisation LU en fonction de la taille
import plotly.graph_objects as go
fig = go.Figure()
fig.add_trace(go.Scatter(x=Ns, y=Ts, mode='markers+lines', name='LU'))
fig.update_layout(xaxis_type="log", yaxis_type="log")
fig.show()
Results

Nous pouvons enfin calculer la pente de la courbe en utilisant la fonction linregress de la bibliothèque scipy.stats

from scipy.stats import linregress
slope, intercept, r_value, p_value, std_err = linregress(np.log(Ns[1:]), np.log(Ts[1:]))
print(f"La pente de la courbe est {slope:.4f}")
Results 
La pente de la courbe est 2.1047