Partiamo dalla definizione che è facile ed intuitiva. Si definisce inversa di una matrice quadrata la matrice che moltiplicata per la matrice stessa restituisce la matrice identità:
Nel caso di matrice di lato unitario è tutto troppo semplice, nel caso di matrice con lato 2 si trova facilmente (si lascia la dimostrazione al volenteroso lettore ) che data:
la matrice inversa è la seguente:
Si nota fin da subito un collegamento con il determinante di una matrice 2x2: le matrici 2x2 per essere invertibili debbono avere determinante non nullo.
La domanda che sorge spontanea è, come faccio a sapere preventivamente se una matrice quadrata di lato maggiore di due è o meno invertibile?
Qui la dimostrazione è decisamente meno facile, si arriva che la matrice deve avere le seguenti proprietà per essere invertibile:
determinante non nullo (determinante diverso da zero);
rango massimo (righe/colonne linearmente indipendenti, vedi sottoparagrafo dedicato al rango);
gli autovalori sono non nulli (autovalori tutti diversi da zero).
Vi sono altre caratteristiche equivalenti a quelle sopra riportate che identificano il gruppo delle matrici invertibili, non è mia intenzione addentrarmi ulteriormente nell'argomento.
Passiamo senza indugio invece all’algoritmo di calcolo classico che è detto metodo di Gauss-Jordan che come spesso accade quando si ha a che fare con Gauss è geniale, il metodo è composto da tre passi sotto riportati.
Prima cosa si prende la matrice quadrata A[nxn] e la si affianca ad una matrice identità I[nxn] ottenendo la matrice aumentata B[nx2n]:
Risulta utile partire subito con un esempio.
Si proceda applicando a B[nxn] il metodo di eliminazione di Gauss che riduce la matrice in una nuova matrice triangolare superiore U[nxn].
Per ricavare U si possono usare solo le seguenti operazioni elementari sulle righe della matrice:
scambiare due righe;
moltiplicare una riga per uno scalare (reale non nullo);
sommare ad una riga un’altra riga moltiplicata per uno scalare.
Per ridurre ad una matrice triangolare superiore la matrice aumentata di partenza si cerca la prima riga con primo elemento non nullo (eventualmente effettuando uno scambio di riga).
Se tutte le righe hanno primo coefficiente nullo si considera la sottomatrice che parte dalla colonna successiva.
Trovato il primo elemento non nullo (detto pivot) tutti gli elementi inferiori della colonna si debbono annullare, per fare ciò basta semplicemente moltiplicare le righe inferiori per un opportuno scalare per la prima riga.
Fatto questo si considera la sotto-matrice che va dalla riga e colonna successive e si reitera il procedimento.
E’ più complicato a dirsi che a fasi, da notare che la matrice risultante non è sempre la stessa e dipende dalle scelte fatte.
Qui si capisce molto meglio con l’esempio.
Prendiamo la prima riga della matrice A, il numero è diverso da 0, il pivot della prima riga quindi risulta:
Procediamo quindi ad annullare gli elementi della prima colonna partendo dalla seconda riga.
Si ha che:
Si passa poi alla terza riga:
E così via per la quarta riga e la quinta riga.
Si otterrà alla fine della computazione la seguente matrice.
Noterete il pivot in alto a destra e sotto nella prima colonna tutti gli zeri.
Ora consideriamo la sotto matrice che si ha eliminando prima riga e prima colonna e reiteriamo il processo, poi prendiamo la sotto matrice ulteriore e avanti così.
Alla fine dopo otteniamo quanto cercato, ovvero una matrice triangolare superiore (nella diagonale ci sono i pivot).
Ora vogliamo far diventare la matrice di cui sopra (triangolare superiore) una matrice così composta:
Si dimostra che la matrice V è proprio l’inversa cercata.
Per ricavare I si procede dall’ultima riga della matrice alla prima e si utilizzano esattamente le stesse tecniche del passo precedente per annullare i valori della colonna sopra il pivot, successivamente moltiplicando per il reciproco dei pivot (che quindi deve essere non nullo) ricaviamo quanto cercato e dunque l’inversa.
Nell'esempio si ricaverà la matrice seguente dove sono ben distinguibili la matrice identità e l'inversa.
Vantaggio non da poco, il metodo è facilmente traducibile in un algoritmo.
Per la risoluzione di inverse di matrici dalle dimensioni generose, tipiche ad esempio delle applicazioni per intelligenza artificiale, esistono algoritmi che stimano l’inversa che sono decisamente più efficienti.
Riporto un brevissimo esempio di utilizzo della classe da me creata per il calcolo dell'inversa della matrice.
//Delcare matrix
Matrix A("[ 1, 0, 3, 0, 5 ; 2, 4, 0, 4, 6 ; 1, 0, 4, 0, 6 ; 1, 0, 2, 4, 3 ; 5, 4, 4, 1, 1 ]");
//Print matrix
std::cout << "A:\n";
print(A);
std::cout << "\n";
//Inverse matrix
Matrix B;
B = A^-1; //Inverse of a square matrix
//Print inverse matrix
std::cout << "A^-1:\n";
print(B);
std::cout << "\n";
//Check A*B
Matrix AB;
AB = A * B;
std::cout << "A * A^-1:\n";
print(AB);
std::cout << "\n";
//Check B*A
Matrix BA;
BA = B * A;
std::cout << "A^-1 * A:\n";
print(BA);
std::cout << "\n";
L'output (con aggiunte anche le matrici intemedie di calcolo viste sopra) sarà il seguente:
A:
1.000 0.000 3.000 0.000 5.000
2.000 4.000 0.000 4.000 6.000
1.000 0.000 4.000 0.000 6.000
1.000 0.000 2.000 4.000 3.000
5.000 4.000 4.000 1.000 1.000
AI=[A|I]:
1.000 0.000 3.000 0.000 5.000 1.000 0.000 0.000 0.000 0.000
2.000 4.000 0.000 4.000 6.000 0.000 1.000 0.000 0.000 0.000
1.000 0.000 4.000 0.000 6.000 0.000 0.000 1.000 0.000 0.000
1.000 0.000 2.000 4.000 3.000 0.000 0.000 0.000 1.000 0.000
5.000 4.000 4.000 1.000 1.000 0.000 0.000 0.000 0.000 1.000
U (gaussian elimination):
1.000 0.000 3.000 0.000 5.000 1.000 0.000 0.000 0.000 0.000
0.000 4.000 -6.000 4.000 -4.000 -2.000 1.000 0.000 0.000 0.000
0.000 0.000 1.000 0.000 1.000 -1.000 0.000 1.000 0.000 0.000
0.000 0.000 0.000 4.000 -1.000 -2.000 0.000 1.000 1.000 0.000
0.000 0.000 0.000 0.000 -15.750 -9.500 -1.000 5.750 0.750 1.000
[I|A^-1]:
1.000 0.000 0.000 0.000 0.000 2.794 -0.127 -2.270 0.095 0.127
0.000 1.000 0.000 0.000 0.000 -1.952 0.202 1.524 -0.214 0.048
0.000 0.000 1.000 0.000 0.000 -1.603 -0.063 1.365 0.048 0.063
0.000 0.000 0.000 1.000 0.000 -0.349 0.016 0.159 0.238 -0.016
0.000 0.000 0.000 0.000 1.000 0.603 0.063 -0.365 -0.048 -0.063
A^-1:
2.794 -0.127 -2.270 0.095 0.127
-1.952 0.202 1.524 -0.214 0.048
-1.603 -0.063 1.365 0.048 0.063
-0.349 0.016 0.159 0.238 -0.016
0.603 0.063 -0.365 -0.048 -0.063
A * A^-1:
1.000 0.000 -0.000 0.000 0.000
-0.000 1.000 0.000 0.000 0.000
0.000 0.000 1.000 0.000 0.000
0.000 0.000 0.000 1.000 0.000
0.000 0.000 0.000 0.000 1.000
A^-1 * A:
1.000 0.000 -0.000 0.000 0.000
-0.000 1.000 0.000 0.000 0.000
0.000 0.000 1.000 0.000 0.000
0.000 0.000 0.000 1.000 0.000
0.000 0.000 0.000 0.000 1.000