Matrice inversa (con decomposizione LU)

Problema: Data una matrice quadrata $\mathbf{A}$  $n \times n$ invertibile, si calcoli la sua inversa $\mathbf{A}^{-1}$.
Soluzione: Sfrutta la decomposizione $\mathbf{LU}$
Complessità: $O(n^3)$

La decomposizione $\mathbf{LU}$ non è altro che una ottimizazione dell'eliminazione di Gauss. Quest'ultima prevede che si operi sulla matrice completa (ottenuta unendo matrice dei coefficienti e vettore dei termini noti) quindi se la matrice dei coefficienti rimane la stessa ma si ha necessità di trovare soluzioni diverse per vari vettori di termini noti si dovrà ripetere tutto il processo di elaborazione ogni volta. Con la decomposizione $\mathbf{LU}$, invece, si opera solo sulla matrice dei coefficienti. In questo caso, anche se si ha necessità di calcolare più soluzioni non si dovrà ripetere nuovamente tutto il processo di elaborazione per ogni nuovo vettore dei termini noti.

Se una matrice $\mathbf{A}$ è invertibile, la si può sempre fattorizzare in questo modo

\[\mathbf{A} = \mathbf{LU}\]

con $\mathbf{L}$ matrice triangolare bassa, che conserverà i fattori moltiplicativi dell'eliminazione di Gauss eseguita su $\mathbf{A}$ (concetto ripreso più avanti) e con $\mathbf{U}$ matrice triangolare alta, che altro non è che il risultato dell'eliminazione di Gauss sulla matrice $\mathbf{A}$.

\[A=\begin{bmatrix} a_{11} & a_{12} & \dots & a_{1n}\\ a_{21} & a_{22} & \dots & a_{2n}\\ \dots  & \dots  & \dots & \dots\\ a_{n1} & a_{n2} & \dots & a_{nn}\\ \end{bmatrix} \quad\quad L=\begin{bmatrix} 1 & 0 & \dots & 0\\ l_{21} & 1 &\dots & 0\\ \dots & \dots & \dots & \dots\\ l_{n1} & l_{n2} & \dots & 1\\ \end{bmatrix} \quad\quad U=\begin{bmatrix} u_{11} & u_{12} & \dots & u_{1n}\\ 0 & u_{22} & \dots & u_{2n}\\ \dots & \dots & \dots & \dots\\ 0 & 0 & \dots & u_{nn}\\ \end{bmatrix}\]

Gli elementi $l_{i,j}$ nella matrice $\mathbf{L}$, come accennato, sono i fattori moltiplicativi usati per modificare le righe della matrice $\mathbf{A}$ quando si effettua l'eliminazione di Gauss. Quindi, ad esempio, se per eliminare l'elemento $a_{2,1}$ nella matrice $\mathbf{A}$ bisogna moltiplicare la prima riga per il fattore $f_{2,1} = a_{2,1}/a_{1,1}$, allora $l_{2,1} = f_{2,1}$.

Abbiamo detto che la decomposizione $\mathbf{LU}$ altro non è che una ottimizzazione dell'eliminazione di Gauss quindi può essere utilizzata anch'essa come tecnica per risolvere sistemi di equazioni lineari che, in forma matriciale, sono del tipo

\begin{equation}\label{eq:1}\mathbf{Ax} = \mathbf{b}\end{equation}

Procedendo con l'eliminazione di Gauss su tale equazione si ottiene

\begin{equation}\label{eq:2}\mathbf{Ux} = \mathbf{d}\end{equation}

dove $\mathbf{U}$ è la matrice triangolare alta ridotta a scala dall'eliminazione di Gauss eseguita su $\mathbf{A}$ e che è proprio il fattore $\mathbf{U}$ della decomposizione $\mathbf{LU}$.
Assumiamo ora che esista una matrice triangolare bassa che se moltiplicata ad entrambi i lati di \eqref{eq:2} da come risultato \eqref{eq:1}. Cioè

\begin{equation}\label{eq:3}\mathbf{LUx} = \mathbf{Ld}\end{equation}

con

\begin{equation}\label{eq:4}\mathbf{LU} = \mathbf{A}\end{equation}

e

\begin{equation}\label{eq:5}\mathbf{Ld} = \mathbf{b}\end{equation}


Guardando la \eqref{eq:4} in effetti tale matrice esiste ed è proprio la matrice triangolare bassa $\mathbf{L}$ i cui elementi sono i fattori moltiplicativi dell'eliminazione di Gauss.
Per risolvere il sistema, dunque, bisogna procedere con un doppio passaggio:

1- Una volta che $\mathbf{L}$ e $\mathbf{b}$ sono noti, risolvere \eqref{eq:5} per trovare $\mathbf{d}$
2- Sostituire in \eqref{eq:2} il vettore trovato al primo passaggio. Ora $\mathbf{U}$ e $\mathbf{d}$ sono noti e si può risolvere \eqref{eq:2} per trovare $\mathbf{x}$

Quindi, per trovare $\mathbf{x}$ nell'equazione \eqref{eq:1} basterà fattorizzare $\mathbf{A}$ con la decomposizione $\mathbf{LU}$ ed eseguire il doppio passaggio visto sopra.
Ed è proprio su questo metodo che si basa il calcolo della matrice inversa. Basterà risolvere

\begin{equation}\label{eq:6}\mathbf{Ax_i} = \mathbf{b_i}\end{equation}

$n$ volte, con $1 \leq i\leq n$. Infatti, impostando $\mathbf{b_i}$ alla $i$-esima colonna di $\mathbf{I}$ (la matrice unitaria) il risultato è che si otterranno in $\mathbf{x_i}$ le colonne della matrice inversa $\mathbf{A^{-1}}$. Questo perché, come è noto, l'unica matrice che moltiplicata ad $\mathbf{A}$ da la matrice unitaria è l'inversa di $\mathbf{A}$.
Le $n$ equazioni da risolvere per trovare le colonne di $\mathbf{A^{-1}}$ si possono sintetizzare in un'unica equazione matriciale

\[\mathbf{AX} = \mathbf{I}\]

Ed è quindi chiaro che $\mathbf{X}$ non può essere altro che $\mathbf{A^{-1}}$. E' proprio questa ripetitività nel risolvere equazioni del tipo \eqref{eq:6} che consente alla decomposizione $\mathbf{LU}$ di trovare l'inversa di una matrice in maniera rapida. Il grosso dell'elaborazione sta proprio nella decomposizione di $\mathbf{A}$. Fatto questo (una sola volta e all'inizio), risolvere le $n$ equazioni del tipo \eqref{eq:6} è relativamente efficiente a livello computazionale.

Prima di vedere come fattorizzare una matrice invertibile con la decomposizione $\mathbf{LU}$ è utile vedere più nel dettaglio i due passaggi esaminati in precedenza per risolvere il sistema di $n$ equazioni. Nel primo passaggio si ha

\[\mathbf{Ld} = \mathbf{b}\]

e cioè

\[\begin{bmatrix} 1 & 0 & \dots & 0\\ l_{21} & 1 & \dots & 0\\ \dots & \dots & \dots & \dots\\ l_{n1} & l_{n2} & \dots & 1\\ \end{bmatrix} \times \mathbf{d} = \mathbf{b}\]

Procedendo con una sostituzione in avanti si possono trovare le componenti del vettore $\mathbf{d}$

\begin{equation}\label{eq:7}d_i = \begin{cases} b_1, & \text{con $i = 1$,}\\ b_i - \sum_{j=1}^{i-1} l_{ij}d_j, &  \text{con $2 \leq i\leq n$.} \end{cases}\end{equation}

Ottenuto $\mathbf{d}$ si può procedere con il secondo passaggio


\[\mathbf{Ux} = \mathbf{d}\]

e cioè

\[\begin{bmatrix} u_{11} & u_{12} & \dots & u_{1n}\\ 0 & u_{22} & \dots & u_{2n}\\ \dots & \dots & \dots & \dots\\ 0 & 0 & \dots & u_{nn}\\ \end{bmatrix}\times \mathbf{x} = \mathbf{d}\]

Procedendo con una sostituzione all'indietro si possono trovare le componenti del vettore $\mathbf{x}$

\begin{equation}\label{eq:8}x_i = \begin{cases} d_n/u_{n,n}, & \text{con $i = n$,}\\ \\ \displaystyle\frac{d_i - \sum_{j=i+1}^{n}  u_{ij}x_j}{u_{ii}}, &  \text{con $1 \leq i < n$.} \end{cases}\end{equation}

Non resta che vedere come implementare in codice la fattorizzazione di una matrice $\mathbf{A}$ attraverso la decomposizione $\mathbf{LU}$. Ricordando cosa rappresentano $\mathbf{U}$ ed $\mathbf{L}$ si può abbozzare un primo tentativo di codice

void lu_decompose(float *lu, int n)
{
 float factor;
 
 // Seleziona k-esima riga
 for (int k = 0; k < n - 1; ++k)
 {
  // Seleziona i-esima riga
  for (int i = k + 1; i < n; ++i)
  {
   // factor = a[i][k]/a[k][k]: fattore moltiplicativo con cui viene moltiplicata la 
   // k-esima riga prima di essere sottratta alla i-esima riga
   factor = lu[n*i + k] / lu[n*k + k];
 
   // Salva fattore nella matrice lu alla posizione (i,k).
   // Poichè i > k, (i,k) corrisponde sempre ad un elemento della parte triangolare bassa della matrice lu
   lu[n*i + k] = factor;
 
   // Eliminazione di Guass per calcolare matrice triangolare alta U.
   // Usa sempre matrice lu dato che non c'è conflitto ad esclusione di diagonale principale (una
   // matrice è triangolare bassa e l'altra è triangolare alta). Poichè i <= j, (i,j) corrisponde 
   // sempre a parte triangolare alta della matrice lu.
   // Indice j seleziona colonna della i-esima riga (o k-esima nel caso della moltiplicazione con factor)
   for (int j = k + 1; j < n; ++j)
    lu[n*i + j] = lu[n*i + j] - factor * lu[n*k + j];
  }
 }
}

Mentre, seguendo quanto visto in \eqref{eq:7} e \eqref{eq:8}, il codice per calcolare i vettori $\mathbf{d}$ e $\mathbf{b}$ può essere scitto in questo modo

void substitute(float *a, int n, float *b, float *x)
{
 float sum;
 
 // Sostituzione in avanti per trovare d; 
 // Segue formula (7) ma con indici che partono da 0 invece che da 1.
 // Usa b come output per d
 // i parte da 1 poichè b[0] resta invariato: d[0] = b[0]
 for (int i = 1; i < n; ++i)
 {
  sum = 0;
 
  // In (7) alla i-esima riga somma (i-1) termini perchè i partiva da 2. 
  // Qui somma i primi i termini perchè i parte da 1
  for (int j = 0; j < i; ++j)
  {
   sum += a[n*i + j] * b[j];
  }
 
  b[i] = b[i] - sum;
 }
 
 // Sostituzione all'indietro per trovare x.
 // Segue formula (8) ma con indici che partono da 0 invece che da 1.
 x[n - 1] = b[n - 1] / a[n*(n - 1) + (n - 1)];  // denom. equivalente a (*(*(a+(n-1))+(n-1)));
 
 for (int i = n - 2; i >= 0; --i)
 {
  sum = 0;
 
  for (int j = i + 1; j < n; ++j)
  {
   sum += a[n*i + j] * x[j];
  }
 
  x[i] = (b[i] - sum) / a[n*i + i];
 }
}

Il codice scritto finora funziona ma adesso bisogna prendere in considerazione alcuni casi che possono portare a risultati non corretti (nel migliore dei casi). Ad esempio, se sulla diagonale principale di $\mathbf{A}$ c'è uno zero, il codice di entrambe le funzioni cercherà di eseguire una divisione con denominatore uguale a 0, con esiti disastrosi. Errori meno catastrofici si verificano a causa della rappresentazione limitata dei numeri reali nei computer. In questo ambito, ad esempio, se si somma un numero a virgola mobile molto grande con uno molto piccolo, vicino allo zero, ci sono grosse probabilità che tale somma non sia rappresentabile sul nostro sistema, causando un errore di arrotondamento. Se all'inizio possono sembrare errori irrilevanti, man mano che i calcoli si fanno più numerosi la situazione si aggrava fino a portare a risultati completamente errati. Per arginare questo tipo di inconvenienti si può ricorrere al pivoting e allo scaling.

Pivoting
Come già accennato, se sulla diagonale principale c'è uno zero si verificherà una eccezione per via della divisione eseguita in entrambe le funzioni. Anche valori vicini a zero sulla diagonale potrebbero dare problemi dato che dividere un numero relativamente grande per un numero molto piccolo produce numeri enormi che potrebbero non essere rappresentabili sulla nostra macchina. Il pivoting può alleviare questi problemi agendo durante l'eliminazione di Gauss scambiando semplicemente le righe della matrice in modo da portare l'elemento più grande di una colonna sulla diagonale principale.
Questo ha l'effetto anche di ridurre gli errori di arrotondamento dato che dividere per numeri grandi produce numeri piccoli ed i floating point hanno una maggiore densità di rappresentazione in questo ambito.
L'implementazione del pivoting prevede scambi tra le righe della matrice ma si ottengono prestazioni migliori anche solo tenendo traccia delle righe da scambiare attraverso un vettore di indici.

Scaling
Durante l'eliminazione di Gauss ci sono varie somme e sottrazioni tra le righe. Se gli elementi delle righe della matrice hanno ordini di grandezza differenti, come già accennato, possono verificarsi problemi di arrotondamento causato dalla limitata rappresentabilità dei numeri reali sui computer. Per evitare somme e sottrazioni tra numeri molto grandi e numeri molto piccoli (tendenti a zero) si scalano tutte le righe della matrice dividendole per l'elemento più grande di ciascuna riga. L'elemento più grande di ogni riga sarà quindi 1 e tutta l'operazione di riduzione della matrice avverrà con numeri relativamente piccoli, dove la rappresentabilità dei floating point ha densità maggiore.
Un effetto collaterale dello scaling, purtroppo, è che manipolare numeri così piccoli attraverso varie operazioni aritmetiche può portare a problemi di rappresentabilità se si va verso valori che tendono quasi allo zero.
Ci si potrebbe domandare a questo punto quale sia l'utilità dello scaling. In realtà scalare le righe è essenziale solo se si vuole calcolare il determinante della matrice per verificare che il sistema non sia mal condizionato. In caso contrario se ne può fare a meno.
Un'implementazione non interessata al calcolo del determinante può in ogni caso usare lo scaling ai soli fini del pivoting, per capire cioè se scambiare o meno due righe senza scalarle realmente.

Tolleranza
Il pivoting, come detto, può aiutare ad alleviare i problemi derivanti da elementi troppo piccoli sulla diagonale principale riordinando le righe. Il problema però persisterebbe se, ad esempio, la colonna presa in esame dal pivoting contiene tutti valori molto vicini allo zero. E' utile quindi usare un valore sentinella durante i calcoli per controllare che sulla diagonale principale ci siano valori accettabili.

Determinante
Calcolare il determinante di una matrice $\mathbf{A}$ è banale quando si ha accesso alla sua decomposizione $\mathbf{LU}$. Il determinante di una matrice triangolare è il prodotto degli elementi della diagonale principale e, per quanto visto fino ad ora, la matrice $\mathbf{U}$ altro non è che la matrice triangolare alta ricavata dalla riduzione di $\mathbf{A}$ con l'eliminazione di Gauss. Basterà quindi moltiplicare gli elementi sulla diagonale principale di $\mathbf{U}$.

Inversa di una matrice
Come già accennato, per trovare l'inversa di una matrice $\mathbf{A}$ basta risolvere $n$ volte l'equazione \eqref{eq:6} per trovare le colonne della sua inversa $\mathbf{A^{-1}}$.

/*Per facilitare la stesura, la lettura e la comprensione del codice si è volutamente 
deciso di utilizzare un mix di C/C++ cercando di restare il più possibile all'interno 
del sottoinsieme C del C++ e di utilizzare solo quei costrutti e funzionalità di base 
del C++ che permettono una maggiore semplificazione.*/
 
#include <stdio.h>
#include <math.h>
 
void inverse(float*, float*, float*, float*, int, int*, float*, float, int&);
void lu_decompose(float*, int, int*, float*, float, int&);
void pivot(float*, int*, float*, int, int);
void substitute(float*, int*, int, float*, float*);
void print_matrix(float*, int);
 
int main(int argc, char *argv[], char *envp[])
{
 float m3x3[3][3] = { {3, -0.1, -0.2}, {0.1, 7, -0.3}, {0.3, -0.2, 10} };
 float tmp3x3[3][3];
 float b[3];
 float x[3];
 int o[3];
 float s[3];
 int err = 0;
 
 // Necessario copiare elementi in matrice di output tmp3x3 perchè verranno modificati solo 
 // gli elementi raggiunti dall'eliminazione di Guass
 for (int i = 0; i < 3; ++i)
 {
  for (int j = 0; j < 3; ++j)
   tmp3x3[i][j] = m3x3[i][j];
 }
 
 inverse(&m3x3[0][0], &tmp3x3[0][0], b, x, 3, o, s, 1e-6, err);
 print_matrix(&tmp3x3[0][0], 3);
 
 return 0;
}
 
void inverse(float *a, float *ai, float *b, float *x, int n, int *o, float *s, float tol, int &err)
{
 // Fattorizza la matrice
 lu_decompose(a, n, o, s, tol, err);
 
 if (!err)
 {
  for (int i = 0; i < n; ++i)
  {
   for (int j = 0; j < n; ++j)
   {
    // Imposta b all'i-esimo vettore della base canonica
    if (i == j)
     b[j] = 1;
    else
     b[j] = 0;
   }
   // Calcola la i-esima colonna della matrice inversa
   substitute(a, o, n, b, x);
 
   // Indici invertiti dato che vettore x contiene colonne di inversa
   for (int j = 0; j < n; ++j)
    ai[n*j + i] = x[j];
  }
 }
 else
  fprintf(stderr, "Errore di tolleranza");
}
 
 
void lu_decompose(float *lu, int n, int *o, float *s, float tol, int &err)
{
 float factor;
 // Inizializza vettore o degli indici delle righe da scambiare e
 // mette nel vettore s l'elemento più grandi di ogni riga
 for (int i = 0; i < n; ++i)
 {
  o[i] = i;
  s[i] = fabs(lu[n*i]); // interessa il valore assoluto
  for (int j = 1; j < n; ++j)
  {
   if (fabs(lu[n*i + j]) > s[i])
    s[i] = fabs(lu[n*i + j]);
  }
 }
 
 // Seleziona k-esima riga attraverso vettore o
 for (int k = 0; k < n - 1; ++k)
 {
  // "Riordina" le righe; in realtà riordina solo un vettore degli indici relativi alle righe
  pivot(lu, o, s, n, k);
  // Controllo di tolleranza sulla diagonale principale.
  // Usa valori scalati senza che le righe vengano scalate davvero.
  if (fabs(lu[n*o[k] + k] / s[o[k]]) < tol)
  {
   err = 1;
   return;
  }
  // Seleziona i-esima riga attraverso o
  for (int i = k + 1; i < n; ++i)
  {
   // factor = a[i][k]/a[k][k]: fattore moltiplicativo con cui viene moltiplicata la 
   // k-esima riga prima di essere sottratta alla i-esima riga
   factor = lu[n*o[i] + k] / lu[n*o[k] + k];
 
   // Salva fattore nella matrice lu alla posizione (i,k).
   // Poichè i > k, (i,k) corrisponde sempre ad un elemento della parte triangolare bassa della matrice lu
   lu[n*o[i] + k] = factor;
 
   // Eliminazione di Guass per calcolare matrice triangolare alta U.
   // Usa sempre matrice lu dato che non c'è conflitto ad esclusione di diagonale principale (una
   // matrice è triangolare bassa e l'altra è triangolare alta). Poichè i <= j, (i,j) corrisponde 
   // sempre a parte triangolare alta della matrice lu.
   // Indice j seleziona colonna della i-esima riga (o k-esima nel caso della moltiplicazione con factor)
   for (int j = k + 1; j < n; ++j)
    lu[n*o[i] + j] = lu[n*o[i] + j] - factor * lu[n*o[k] + j];
  }
 }
 // Controllo di tolleranza su ultimo elemento di diagonale principale.
 // Usa valori scalati senza che la riga venga scalata davvero
 if (fabs(lu[n*o[n - 1] + (n - 1)] / s[o[n - 1]]) < tol)
  err = 1;
}
 
 
void pivot(float *a, int *o, float *s, int n, int k)
{
 // p manterrà l'indice con l'elemento più grande all'interno della stessa colonna.
 // big manterrà l'elemento più grande all'interno della stessa colonna.
 // Usa valori scalati per valutare lo scambio ma senza scalare le righe per davvero.
 int p = k;
 float big = fabs(a[n*o[k] + k] / s[o[k]]);
 float dummy;
 
 // Seleziona i-esima riga nella k-esima colonna
 for (int i = k + 1; i < n; ++i)
 {
  dummy = fabs(a[n*i + k] / s[o[i]]);
  if (dummy > big)
  {
   big = dummy;
   p = i;
  }
 }
 // "Scambia" p-esima riga con k-esima riga.
 // In realtà scambia solo gli indici all'interno del vettore o
 int z = o[p];
 o[p] = o[k];
 o[k] = z;
}
 
 
void substitute(float *a, int *o, int n, float *b, float *x)
{
 float sum;
 
 // Sostituzione in avanti per trovare d; 
 // Segue formula (7) ma con indici che partono da 0 invece che da 1.
 // Usa b come output per d
 // i parte da 1 poichè b[0] resta invariato: d[0] = b[0]
 
 // In (7) alla i-esima riga somma (i-1) termini perchè i partiva da 2. 
 // Qui somma i primi i termini perchè i parte da 1
 for (int i = 1; i < n; ++i)
 {
  sum = 0;
  for (int j = 0; j < i; ++j)
  {
   sum += a[n*o[i] + j] * b[o[j]];
  }
  b[o[i]] = b[o[i]] - sum;
 }
 
 // Sostituzione all'indietro per trovare x.
 // Segue formula (8) ma con indici che partono da 0 invece che da 1.
 x[n - 1] = b[o[n - 1]] / a[n*o[(n - 1)] + (n - 1)];
 for (int i = n - 2; i >= 0; --i)
 {
  sum = 0;
  for (int j = i + 1; j < n; ++j)
  {
   sum += a[n*o[i] + j] * x[j];
  }
  x[i] = (b[o[i]] - sum) / a[n*o[i] + i];
 }
}
 
void print_matrix(float *m, int n)
{
 for (int i = 0; i < n; ++i)
 {
  for (int j = 0; j < n; ++j)
  {
   printf("%f ", m[n*i + j]);
  }
 
  printf("\n");
 }
 printf("\n");
}

Infine, questo è il confronto tra il risultato del nostro codice e quello restituito da Wolfram|Alpha

Riferimenti
[1] Numerical Methods for Engineers - Chapra, Canale
[2] Numerical Recipies - Press, Teukolsky, Vetterling, Flannery
[3] Introduction to Algorithms - Cormen, Leiserson, Rivest, Stein