Eliminarea Gauss-Jordan#

Asemănător eliminării gaussiane, mai există un algoritm care prezintă interes inginerilor, denumit Eliminare Gauss-Jordan.

Acesta are două versiuni, una care îndeplinește scopul obișnuit al eliminării gaussiene, anume soluționarea unui sistem de ecuații liniare, respectiv una care este utilizată pentru a inversa o matrice.

Ambele vor fi discutate amănunțit în acest subcapitol.

Vizualizarea algoritmului#

Spre deosebire de algoritmii clasici de eliminare Gaussiană, algoritmul de eliminare Gauss-Jordan își propune transformarea matricei \(A\in\mathbb{R}^{n\times n}\), \(n\in\mathbb{N}^*\), într-o matrice identitate (\(I_n\)), nu doar superior triunghiulară.

Să exemplificăm pentru o matrice de coeficienți \(A\in\mathbb{R}^{4\times4}\) ce descrie sistemul de ecuații liniare \(A\bm{x}=\bm{b}\), unde \(\bm{x},\bm{b}\in\mathbb{R}^4\):

\[ A\bm{x}=\bm{b}\Leftrightarrow\begin{bmatrix} \boxtimes & \boxtimes & \boxtimes & \boxtimes \\ \boxtimes & \boxtimes & \boxtimes & \boxtimes \\ \boxtimes & \boxtimes & \boxtimes & \boxtimes \\ \boxtimes & \boxtimes & \boxtimes & \boxtimes \end{bmatrix} \begin{bmatrix} x_1\\x_2\\x_3\\x_4 \end{bmatrix}=\begin{bmatrix} b_1\\b_2\\b_3\\b_4 \end{bmatrix} \]

Vrem ca elementele de pe diagonala principală să fie egale cu \(1\), așadar împărțim prima linie la valoarea \(a_{11}\) și obținem:

\[ \begin{bmatrix} a_{11} & \boxtimes & \boxtimes & \boxtimes \\ \boxtimes & \boxtimes & \boxtimes & \boxtimes \\ \boxtimes & \boxtimes & \boxtimes & \boxtimes \\ \boxtimes & \boxtimes & \boxtimes & \boxtimes \end{bmatrix} \begin{bmatrix} x_1\\x_2\\x_3\\x_4 \end{bmatrix}=\begin{bmatrix} b_1\\b_2\\b_3\\b_4 \end{bmatrix}\rightarrow\begin{bmatrix} 1 & \boxtimes & \boxtimes & \boxtimes \\ \boxtimes & \boxtimes & \boxtimes & \boxtimes \\ \boxtimes & \boxtimes & \boxtimes & \boxtimes \\ \boxtimes & \boxtimes & \boxtimes & \boxtimes \end{bmatrix} \begin{bmatrix} x_1\\x_2\\x_3\\x_4 \end{bmatrix}=\begin{bmatrix} b_1'\\b_2\\b_3\\b_4 \end{bmatrix} \]

Vom elimina toate elementele care se află pe prima coloană și nu se regăsesc pe diagonala principală, așa cum a fost descris procesul în cadrul eliminării gaussiene clasice:

\[ \begin{bmatrix} 1 & \boxtimes & \boxtimes & \boxtimes \\ \boxtimes & \boxtimes & \boxtimes & \boxtimes \\ \boxtimes & \boxtimes & \boxtimes & \boxtimes \\ \boxtimes & \boxtimes & \boxtimes & \boxtimes \end{bmatrix} \begin{bmatrix} x_1\\x_2\\x_3\\x_4 \end{bmatrix}=\begin{bmatrix} b_1'\\b_2\\b_3\\b_4 \end{bmatrix}\rightarrow\begin{bmatrix} 1 & \boxtimes & \boxtimes & \boxtimes \\ 0 & \boxtimes & \boxtimes & \boxtimes \\ 0 & \boxtimes & \boxtimes & \boxtimes \\ 0 & \boxtimes & \boxtimes & \boxtimes \end{bmatrix} \begin{bmatrix} x_1\\x_2\\x_3\\x_4 \end{bmatrix}=\begin{bmatrix} b_1'\\b_2'\\b_3'\\b_4' \end{bmatrix} \]

Procesul poate părea plictisitor până la acest moment, însă insistăm să ilustrăm și cea de-a doua coloană - pentru început, se împarte a doua linie la \(a_{22}'\):

\[ \begin{bmatrix} 1 & \boxtimes & \boxtimes & \boxtimes \\ 0 & a_{22}' & \boxtimes & \boxtimes \\ 0 & \boxtimes & \boxtimes & \boxtimes \\ 0 & \boxtimes & \boxtimes & \boxtimes \end{bmatrix} \begin{bmatrix} x_1\\x_2\\x_3\\x_4 \end{bmatrix}=\begin{bmatrix} b_1'\\b_2'\\b_3'\\b_4' \end{bmatrix}\rightarrow\begin{bmatrix} 1 & \boxtimes & \boxtimes & \boxtimes \\ 0 & 1 & \boxtimes & \boxtimes \\ 0 & \boxtimes & \boxtimes & \boxtimes \\ 0 & \boxtimes & \boxtimes & \boxtimes \end{bmatrix} \begin{bmatrix} x_1\\x_2\\x_3\\x_4 \end{bmatrix}=\begin{bmatrix} b_1'\\b_2''\\b_3'\\b_4' \end{bmatrix} \]

Mergem pe aceeași idee de a elimina elementele ce nu se află pe diagonala principală din cea de-a doua coloană (de data aceasta avem nevoie de mai mult decât un simplu algoritm G) și obținem:

\[ \begin{bmatrix} 1 & \boxtimes & \boxtimes & \boxtimes \\ 0 & 1 & \boxtimes & \boxtimes \\ 0 & \boxtimes & \boxtimes & \boxtimes \\ 0 & \boxtimes & \boxtimes & \boxtimes \end{bmatrix} \begin{bmatrix} x_1\\x_2\\x_3\\x_4 \end{bmatrix}=\begin{bmatrix} b_1'\\b_2''\\b_3'\\b_4' \end{bmatrix}\rightarrow\begin{bmatrix} 1 & 0 & \boxtimes & \boxtimes \\ 0 & 1 & \boxtimes & \boxtimes \\ 0 & 0 & \boxtimes & \boxtimes \\ 0 & 0 & \boxtimes & \boxtimes \end{bmatrix} \begin{bmatrix} x_1\\x_2\\x_3\\x_4 \end{bmatrix}=\begin{bmatrix} b_1''\\b_2''\\b_3''\\b_4'' \end{bmatrix} \]

Observăm deci cum algoritmul acesta nu modifică doar valorile aflate sub diagonala principală, ci modifică întreaga matrice \(A\). Rezultatul final va fi dat de relația:

\[ A\bm{x}=\bm{b}\Leftrightarrow\begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} x_1\\x_2\\x_3\\x_4 \end{bmatrix}=\begin{bmatrix} b_1^{(4)}\\b_2^{(4)}\\b_3^{(4)}\\b_4^{(4)} \end{bmatrix}\Rightarrow\begin{bmatrix} x_1\\x_2\\x_3\\x_4 \end{bmatrix}=\begin{bmatrix} b_1^{(4)}\\b_2^{(4)}\\b_3^{(4)}\\b_4^{(4)} \end{bmatrix} \]

Vom considera această ilustrație validă pentru ambele utilizări:

Algoritmul nemodificat - utilizat pentru a soluționa un sistem de ecuații liniare;
Algoritmul modificat - utilizat pentru a inversa o matrice pătratică.

Algoritmul nemodificat#

Să încercăm să formalizăm vizualizarea anterioară.

Fie \(A\bm{x}=\bm{b}\) un sistem de ecuații liniare, unde \(A\in\mathbb{R}^{n\times n}\) este o matrice pătratică și \(\bm{x},\bm{b}\in\mathbb{R}^n\) sunt vectori coloană. Notăm \(\overline{A}=\begin{bmatrix}A&\bm{b}\end{bmatrix}\) matricea extinsă ce conține elemente de forma \(\overline{a}_{ij}\) pe linia \(i\in\overline{1,n}\), coloana \(j\in\overline{1,n+1}\).

Algoritmul de eliminare Gauss-Jordan poate fi descris iterativ astfel: pentru fiecare pas \(i\in\overline{1,n}\), se interschimbă linia \(i\) a matricei \(\overline{A}\) cu oricare linie \(k\in\overline{i,n}\) care conține un element nenul pe coloana \(i\).

Valoarea aleasă

Este de preferat să se foloseasă valoarea maximă posibilă, însă acest lucru nu este obligatoriu, deși asigură o stabilitate numerică mult mai sporită decât varianta obișnuită.

Apoi, în urma interschimbării, elementul \(\overline{a}_{ii}\) devine \(1\) prin împărțirea întregii linii la \(\overline{a}_{ii}\).

În continuare toate elementele matricei \(\overline{A}\) aflate pe coloana \(i\), cu excepția celui aflat pe diagonala principală, se transformă în \(0\) prin aplicarea unui algoritm tip G.

Matematic, pasul \(i\) arată astfel:

Interschimbarea celor poate fi văzută drept drept:
\[ \exists\,k\in\overline{i,n}:\, \overline{a}_{ki}\neq0\Rightarrow\boxed{ \bm{L_i}(\overline{A})\leftrightarrow\bm{L_k}(\overline{A})} \]
Apoi, elementul \(\overline{a}_{ii}\) devine \(1\) prin împărțirea întregii linii la \(\overline{a}_{ii}\):
\[ \boxed{\bm{L_i}(\overline{A})\rightarrow\frac{1}{\overline{a}_{ii}}\cdot\bm{L_i}(\overline{A})} \]
Apoi, se poate aplica o adaptare a algoritmului G clasic:
\[ \boxed{\bm{L_j}(\overline{A})\rightarrow \bm{L_j}(\overline{A})-\overline{a}_{ji}\cdot\bm{L_i}(\overline{A})}\,,\,\,\forall j\neq i \]

Exemplificarea algoritmului#

Vom reconsidera sistemul de ecuații liniare, descris direct în forma matriceală, \(\overline{A}=\begin{bmatrix} 1 & -2 & -3 & | & 0\\ 3 & 5 & 2 & | & 0\\ 2 & 3 & -1 & | & -1 \end{bmatrix}\).

Urmăm pașii descriși anterior:

Pasul 1a. Prima linie rămâne intactă, întrucât conține pe prima poziție o valoare nenulă;
Pasul 1b. Din nou, prima linie rămâne neafectată, deoarece ar fi necesară o împărțire la \(1\) care nu schimbă nicio valoare a liniei;
Pasul 1c. Se golește prima coloană. Se execută:
\[ \begin{cases} \bm{L_2}(\overline{A})\rightarrow\bm{L_2}(\overline{A})-3\bm{L_1}(\overline{A})\\ \bm{L_3}(\overline{A})\rightarrow\bm{L_3}(\overline{A})-2\bm{L_1}(\overline{A}) \end{cases}\Rightarrow\overline{A}=\begin{bmatrix} 1 & -2 & -3&|&0\\0 & 11 & 11&|&0\\0 & 7 & 5&|&-1 \end{bmatrix} \]
Pasul 2a. A doua linie rămâne nemodificată, deoarece conține pe a doua poziție o valoare nenulă, mai cu seamă cea de modul maxim;
Pasul 2b. Se împarte a doua linie la \(\overline{a}_{22}=11\) - se execută:
\[ \bm{L_2}(\overline{A})\rightarrow\frac{1}{11}\bm{L_2}(\overline{A})\Rightarrow\overline{A}=\begin{bmatrix} 1 & -2 & -3&|&0\\0 & 1 & 1&|&0\\0 & 7 & 5&|&-1 \end{bmatrix} \]
Pasul 2c. Se golește a doua coloană. Se execută:
\[ \begin{cases} \bm{L_1}(\overline{A})\rightarrow\bm{L_1}(\overline{A})+2\bm{L_2}(\overline{A})\\ \bm{L_3}(\overline{A})\rightarrow\bm{L_3}(\overline{A})-7\bm{L_2}(\overline{A}) \end{cases}\Rightarrow\overline{A}=\begin{bmatrix} 1 & 0 & -1&|&0\\0 & 1 & 1&|&0\\0 & 0 & -2&|&-1 \end{bmatrix} \]
Pasul 3a. Ultima linie rămâne intactă de asemenea, deoarece conține o valoare nenulă, ba chiar de modul maxim;
Pasul 3b. Se împarte ultima linie la \(\overline{a}_{33}=-2\) - se exceută:
\[ \bm{L_3}(\overline{A})\rightarrow-\frac{1}{2}\bm{L_3}(\overline{A})\Rightarrow\overline{A}=\begin{bmatrix} 1 & 0 & -1&|&0\\0 & 1 & 1&|&0\\0 & 0 & 1&|&0.5 \end{bmatrix} \]

Pasul 3c. Se golește a treia coloană. Se execută:
\[ \begin{cases} \bm{L_1}(\overline{A})\rightarrow\bm{L_1}(\overline{A})-\bm{L_3}(\overline{A})\\ \bm{L_2}(\overline{A})\rightarrow\bm{L_2}(\overline{A})+\bm{L_3}(\overline{A}) \end{cases}\Rightarrow\overline{A}=\begin{bmatrix} 1 & 0 & 0&|&0.5\\0 & 1 & 0&|&-0.5\\0 & 0 & 1&|&0.5 \end{bmatrix} \]

Concluzie. Rezultatul este vectorul coloană din dreapta, adică \(x=\begin{bmatrix} 0.5\\-0.5\\0.5 \end{bmatrix}\).

Cod ilustrativ#

Codul de mai jos implementează algoritmul Gauss-Jordan nemodificat (sub forma funcției Gauss_Jordan), primind parametrii reprezentativi sistemului \(A\bm{x}=\bm{b}\):

A - matricea coeficienților sistemului;
b - vectorul coloană din partea dreaptă a egalității.

Funcția returnează perechea [A, b], echivalentul unui sistem triunghiular de ecuații liniare pe care se poate aplica apoi SST-ul explicat anterior.

function b = Gauss_Jordan(A, b)
    n = length(A);
    Ae = [A, b];
    for p = [1 : n]
        % Folosim GPP pentru ca e mai stabil decat alegerea unui
        % element nenul oarecare:
        [piv, l_piv] = max(abs(Ae(p:n, p)));
        l_piv = l_piv + p - 1;
        temp = Ae(l_piv, :);
        Ae(l_piv,:) = Ae(p, :);
        Ae(p, :) = temp;
        if (Ae(p, p) == 0)
            b = inf;
            return;
        end
        Ae(p, :) /= Ae(p, p);
        for l = 1 : n
            if l != p
                Ae(l, :) -= Ae(l, p) * Ae(p, :);
            end
        end
    end
    b = Ae(:, n+1);
end

Analiza complexității#

Vom considera \(GJ_i(n,m)\) complexitatea algoritmului la fiecare pas \(i\in\overline{1,n}\). Pentru a căuta valoarea lui \(k\) (fie ea aleasă astfel încât să indice un element de modul maxim sau doar o valoare nenulă), sunt necesare cel mult \(n\) comparări.

Folosind un tabel de mapare, interschimbarea a două linii poate fi realizată în timp constant (\(O(1)\)); cu toate acestea, operația costisitoare de a scădea din fiecare element o valoare scalată presupune \((n-1)\cdot m\) operații. Așadar, \(GJ_i(n,m)\in O(mn)\).

Complexitatea finală poate fi obținutăm însumând toate complexitățile \(GJ_i(n,m)\), adică:

\[ GJ(n,m) = \sum_{i=1}^{n}GJ_i(n,m)\in O(mn^2) \]

Algoritmul modificat#

Motivul pentru care algoritmul Gauss-Jordan este impresionant este această modificare a sa care permite inversarea unei matrice pătratice.

Fie \(A\in\mathbb{R}^{n\times n}\), \(n\in\mathbb{N}^*\), o matrice pătratică nesingulară pe care vrem să o inversăm, adică vrem să calculăm \(A^{-1}\). Dacă notăm \(\overline{A}=\begin{bmatrix} A&I_n \end{bmatrix}\), putem refolosi întregul algoritm descris anterior pentru a compune în partea dreaptă matricea \(A^{-1}\).

Demonstrație

Privind matricea \(A\) drept o matrice de coeficienți a unui sistem de forma \(A\bm{x}=\bm{b}\), unde \(\bm{x},\bm{b}\in\mathbb{R}^n\), cunoaștem deja că rezultatul sistemului nu se modifică în urma executării unor transformări elementare asupra lui \(A\).

Așadar, există transformările \(E_1,\dots,E_\mu\), \(\mu\in\mathbb{N}\), alese astfel încât sistemul să se reducă la unul reprezentat de matricea identitate:

\[ E_\mu\cdot E_{\mu-1}\cdot\dots\cdot E_2\cdot E_1\cdot A=I_n \]

Știm însă că, prin aplicarea eliminării Gauss-Jordan, se reduce matricea \(A\) la matricea identitate, așadar cunoaștem algoritmul prin care să aflăm \(E_1,\dots,E_\mu\). Înmulțind cu \(A^{-1}\) în ambele părți, găsim următoarea formulă pentru inversa lui \(A\):

\[ E_\mu\cdot E_{\mu-1}\cdot\dots\cdot E_2\cdot E_1\cdot I_n=A^{-1} \]

În concluzie, urmând algoritmul deja cunoscut până la reducerea matricei \(A\) în matrice identitate se va obține inversa matricei \(A\).

Exemplificarea algoritmului#

Să inversăm matricea \(A=\begin{bmatrix} 2 & 5 & 10\\-1 & -2 & -3\\5 & 6 & 0 \end{bmatrix}\). Pentru aceasta, formăm matricea \(\overline{A}=\begin{bmatrix} 2 & 5 & 10&|&1&0&0\\-1 & -2 & -3&|&0&1&0\\5 & 6 & 0&|&0&0&1 \end{bmatrix}\), iar apoi urmăm algoritmul cunoscut:

Pasul 1a. Prima linie rămâne intactă, întrucât conține pe prima poziție o valoare nenulă;
Pasul 1b. Se împarte prima linie la \(\overline{a}_{11}=2\) - se exceută:
\[ \bm{L_1}(\overline{A})\rightarrow\frac{1}{2}\bm{L_1}(\overline{A})\Rightarrow\overline{A}=\begin{bmatrix} 1 & 2.5 & 5&|&0.5&0&0\\-1 & -2 & -3&|&0&1&0\\5 & 6 & 0&|&0&0&1 \end{bmatrix} \]
Pasul 1c. Se golește prima coloană. Se execută:
\[ \begin{cases} \bm{L_2}(\overline{A})\rightarrow\bm{L_2}(\overline{A})+\bm{L_1}(\overline{A})\\ \bm{L_3}(\overline{A})\rightarrow\bm{L_3}(\overline{A})-5\bm{L_1}(\overline{A}) \end{cases}\Rightarrow\overline{A}=\begin{bmatrix} 1 & 2.5 & 5&|&0.5&0&0\\0 & 0.5 & 2&|&0.5&1&0\\0 & -6.5 & -25&|&-2.5&0&1 \end{bmatrix} \]
Pasul 2a. A doua linie rămâne nemodificată, deoarece conține pe a doua poziție o valoare nenulă, chiar dacă nu cea de modul maxim;
Pasul 2b. Se împarte a doua linie la \(\overline{a}_{22}=0.5\) - se exceută:
\[ \bm{L_2}(\overline{A})\rightarrow2\bm{L_2}(\overline{A})\Rightarrow\overline{A}=\begin{bmatrix} 1 & 2.5 & 5&|&0.5&0&0\\ 0 & 1 & 4&|&1&2&0\\ 0 & -6.5 & -25&|&-2.5&0&1 \end{bmatrix} \]
Pasul 2c. Se golește a doua coloană. Se execută:
\[ \begin{cases} \bm{L_1}(\overline{A})\rightarrow\bm{L_1}(\overline{A})-2.5\bm{L_k}(\overline{A})\\ \bm{L_3}(\overline{A})\rightarrow\bm{L_3}(\overline{A})+6.5\bm{L_k}(\overline{A}) \end{cases}\Rightarrow\overline{A}=\begin{bmatrix} 1 & 0 & -5&|&-2&-5&0\\ 0 & 1 & 4&|&1&2&0\\ 0 & 0 & 1&|&4&13&1 \end{bmatrix} \]
Pasul 3a. Ultima linie rămâne intactă de asemenea, deoarece conține o valoare nenulă, ba chiar de modul maxim;
Pasul 3b. Totodată, a treia linie rămâne neafectată, întrucât ar trebui efectuată o împărțire la \(\overline{a}_{33}=1\) care nu schimbă nicio valoare a liniei.
Pasul 3c. Se golește a treia coloană. Se execută:
\[ \begin{cases} \bm{L_1}(\overline{A})\rightarrow\bm{L_1}(\overline{A})+5\bm{L_3}(\overline{A})\\ \bm{L_k}(\overline{A})\rightarrow\bm{L_2}(\overline{A})-4\bm{L_3}(\overline{A}) \end{cases}\Rightarrow\overline{A}=\begin{bmatrix} 1 & 0 & 0&|&18&60&5\\ 0 & 1 & 0&|&-15&-50&-4\\ 0 & 0 & 1&|&4&13&1 \end{bmatrix} \]

Concluzie. S-a obținut deci \(A^{-1}=\begin{bmatrix} 18 & 60 & 5\\ -15 & -50 & -4\\ 4 & 13 & 1 \end{bmatrix}\).

Cod ilustrativ#

Pentru a implementa varianta modificată a algoritmului Gauss-Jordan prin funcția Gauss_Jordan_inv, vom utiliza notația A pentru matricea pe care încercăm să o inversăm.

Funcția returnează matricea invA, echivalentul lui \(A^{-1}\).

function invA = Gauss_Jordan_inv(A)
    n = length(A);
    Ae = [A, eye(n)];
    for p = [1 : n]
        % Folosim GPP pentru ca e mai stabil decat alegerea unui
        % element nenul oarecare:
        [piv, l_piv] = max(abs(Ae(p:n, p)));
        l_piv = l_piv + p - 1;
        temp = Ae(l_piv, :);
        Ae(l_piv,:) = Ae(p, :);
        Ae(p, :) = temp;
        if (Ae(p, p) == 0)
            b = inf;
            return;
        end
        Ae(p, :) /= Ae(p, p);
        for l = 1 : n
            if l != p
                Ae(l, :) -= Ae(l, p) * Ae(p, :);
            end
        end
    end
    invA = Ae(:, n+1:2*n);
end

Licență#

The book "Metode Numerice", written by Valentin-Ioan Vintilă, is licensed under CC BY-NC-SA 4.0