Matrice Pseudo-inversa

Uno dei problemi fondamentali nell'analisi di sistemi reali è quello del risolvere sistemi lineari sovradimensionati affetti da rumore.

L'importanza di questo argomento é evidente: quando si eseguono osservazioni di un sistema reale questo risulta normalmente affetto da rumore di osservazione (ed eventualmente di processo) ma allo stesso tempo è possibile raccogliere molti più dati che incognite ottenendo normalmente un sistema sovradimensionato. In questa situazione per ottenere una soluzione del problema e allo stesso tempo minimizzare l'errore è richiesto l'utilizzo di una regressione numerica, per esempio ai minimi quadrati.

Si abbia pertanto un sistema lineare sovradimensionato (overdetermined)

$$\mathbf{A} \mathbf{x}= \mathbf{y}$$ (1.1)

dove $\mathbf{A}$ è una matrice rettangolare $m \times n$ e con $m \geq n$. Tale matrice, essendo rettangolare, non ammette inversa ma è comunque possibile definire, per ogni possibile soluzione $\mathbf{x} \in \mathbb{R}^{n}$, un valore dell'errore (residuo) che questa soluzione comporterebbe. Non c'è una soluzione generale per un sistema sovradimensionato, ma solo soluzioni che minimizzano il residuo sotto una particolare norma.

Definiamo, per esempio, come metrica di errore il modulo del residuo:

$$\epsilon(\mathbf{x}) = \left\Vert \mathbf{A} \mathbf{x} - \mathbf{y} \right \Vert^{2}$$ (1.2)

La soluzione ai minimi quadrati di un sistema lineare è rappresenta dal vettore che minimizza la distanza euclidea del residuo 1.2.

Trovare la soluzione ottima del sistema 1.1 nei sensi di una regresione ai minimi quadrati equivale a trovare il minimo di tale funzione errore al variare di $\mathbf{x}$. Va subito fatto notare che nel minimizzare 1.2 non si fa comunque nessuna ipotesi della distribuzione del rumore all'interno dei componenti. La soluzione ottenuta con questa minimizzazione è una soluzione algebrica che minimizza un errore algebrico (algebraic error). Nel capitolo di statistica si affronterà il caso in cui si conosce come l'errore entra nel sistema.

Si può dimostare che una soluzione $\mathbf{x}$, che minimizza la funzione 1.2, esiste e vale:

$$\begin{array}{rcl} \mathbf{A} \mathbf{x} & = & \mathbf{y} ... ...athbf{A}\right)^{-1}\mathbf{A}^{\top} \mathbf{y} \end{array}$$ (1.3)

Per costruzione $\mathbf{x}$ è una soluzione del sistema 1.1 ed è anche il vettore che minimizza la funzione 1.2. Viene indicata con $\mathbf{A}^{+}$ la matrice pseudoinversa (pseudoinverse matrix) di $\mathbf{A}$ e vale

$$\mathbf{A}^{+} = \left(\mathbf{A}^{\top}\mathbf{A}\right)^{-1}\mathbf{A}^{\top}$$ (1.4)

Questa soluzione del sistema è detta pseudoinversa di Moore-Penrose.

La pseudoinversa ha le seguenti proprietà

• La pseudoinversa di una matrice esiste se esiste l'inversa di $\mathbf{A}^{\top}\mathbf{A}$;
• La pseudoinversa di una matrice quadrata coincide con la sua inversa;
• La pseudoinversa di una matrice, se esiste, è unica.

È possibile cercare di migliorare la soluzione in presenza di rumore conosciuto in quanto è possibile assegnare alle equazioni del sistema pesi diversi e questi possono essere rappresentati in una matrice diagonale di precondizionamento della matrice $\mathbf{A}$. Allo stesso modo visto che si minimizza la distanza euclidea della soluzione è possibile moltiplicare ogni riga del sistema per un opportuno peso in modo da pesare in maniera differente ogni dato acquisito.

In lettura è possibile trovare la risoluzione di sistemi sovradimensionati indicata anche come tecnica delle equazioni perpendicolari (normal equations):

$$\mathbf{A}^{\top}\mathbf{A} \mathbf{x} = \mathbf{A}^{\top} \mathbf{y}$$ (1.5)

È facile notare che tale forma è la stessa rappresentazione del sistema 1.3 dove però non viene costruita la matrice pseudoinversa e il problema viene ricondotto a un sistema lineare classico dove la matrice dei coefficienti è quadrata (e pertanto invertibile con tecniche classiche).

La soluzione proposta in equazione 1.5 tuttavia è numericamente instabile in quanto cond$\approx \mathbf{A}^2$. Dettagli ulteriori sul condizionamento delle matrici e sulla propagazione dei disturbi nella soluzione dei sistemi lineari ben dimensionati o sovradimensionati saranno presentati in sezione 2.7.

Se il sistema è ben condizionato, la tecnica più stabile per risolvere un problema alle normal equations è la fattorizzazione di Cholesky.

Esistono invece delle tecniche stabili che permettono di ricavare la soluzione partendo direttamente dalla matrice $\mathbf{A}$. La matrice pseudoinversa si può ottenere, oltre che eseguendo il procedimento mostrato nell'equazione 1.4, anche utilizzando la Singular Value Decomposition (SVD) o la fattorizzazione QR.

Per quanto riguarda la fattorizzazione QR della matrice $\mathbf{A}$ il problema originale 1.1 si trasforma in $\mathbf{Q}\mathbf{R} \mathbf{x} = \mathbf{y}$ e la soluzione si può ricavare da $\mathbf{R}\mathbf{x} = \mathbf{Q}^{\top} \mathbf{y}$.

Attraverso invece la decomposizione ai valori singolari, la matrice sovradimensionata $\mathbf{A}$ viene scomposta in diverse altre matrici con proprietà molto interessanti. Sia $\mathbf{A}=\mathbf{U} \mathbf{S} \mathbf{V}^{*}$ la decomposizione ai valori singolari (SVD) di $\mathbf{A}$. $\mathbf{U}$ è una matrice unitaria di dimensioni $m \times n$ (a seconda del formalismo usato, complete SVD o economic SVD, le dimensioni delle matrici possono cambiare, e $\mathbf{U}$ diventare $m \times m$), $\mathbf{S}$ è una matrice diagonale che contiene i valori singolari (gli autovalori della matrice psudoinversa di dimensioni, a seconda del formalismo , $n \times n$ o $m \times n$) e $\mathbf{V}^{*}$ è una matrice ortonormale, trasposta coniugata, di dimensioni $n \times n$.

Attraverso un procedimento puramente matematico si ottiene che la pseudoinversa di $\mathbf{A}$ equivale a $\mathbf{A}^{+}=\mathbf{V} \mathbf{S}^{+} \mathbf{U}^{*}$ dove la pseudoinversa di una matrice diagonale $\mathbf{S}^{+}$ equivale alla sua inversa ovvero una matrice diagonale costituita dai reciproci dei rispettivi valori.

Dalla SVD si ottiene immediatamente la soluzione

$$\mathbf{x} = \mathbf{A}^{+} \mathbf{y} = \mathbf{V} \mathbf{S}^{-1} \mathbf{U} \mathbf{y}$$ (1.6)

Attraverso la decomposizione SVD è anche possibile trovare lo spazio delle soluzioni di un sistema omogeneo sovradimensionato. Un sistema lineare omogeneo ha la forma

$$\mathbf{A}\mathbf{x}=0$$ (1.7)

e normalmente la soluzione ovvia, anche ottenuta attraverso l'equazione 1.6, non risulta utile ai fini del problema. In questo caso è necessario trovare, sempre ai sensi di una regressione ai minimi quadrati, un $\mathbf{x} \in \mathbb{R}^{n}$, non nullo, ma con un vincolo aggiuntivo per esempio $\vert \mathbf{x} \vert = 1$, tale che

$$\min_{\mathbf{x}} \Vert \mathbf{A}\mathbf{x} \Vert^{2}$$ (1.8)

Anche in questo caso la SVD si dimostra una tecnica estremamente efficiente e computazionalmente stabile, siccome le basi del kernel di $\mathbf{A}$ sono esattamente le colonne di $\mathbf{V}$ associate ai valori (autovalori) nulli della matrice diagonale $\mathbf{S}$. In genere, a causa della presenza di rumore, non esisterà un valore singolare nullo ma si prende solitamente la colonna associata al minimo valore singolare.

Gli autovalori nulli della matrice $\mathbf{S}$ rappresentano pertanto il kernel della matrice stessa. Il numero di autovalori nulli rappresenta la dimensione del kernel stesso. Va notato come nell'equazione 1.6 la presenza di zeri nella matrice $\mathbf{S}$ fosse problematica ma tale presenza è comunque sintomo del fatto che una delle componenti del problema è totalmente uncorrelata con la soluzione e, in quanto tale, può essere trascurata (tale risultato sarà utilizzato in seguito nella trattazione della PCA).

La decomposizione SVD risulta una delle tecniche più stabili e versatili sviluppata negli ultimi anni per la risoluzione di sistemi lineari e di fatto in tutto questo libro si farà larghissimo uso di tale tecnologia.

Dettagli ulteriori sulla pseudoinversa di Moore-Penrose possono essere trovati in molti libri, per esempio in [5] o nel testo fondamentale di calcolo numerico [12].