Sistemi lineari sovradimensionati

Analizzando sistemi reali è facile imbattersi nel problema di dover ricavare la “soluzione” di un sistema lineare sovradimensionato.

L'importanza di questo argomento è evidente: quando si eseguono osservazioni su un sistema reale questo risulta naturalmente affetto da rumore, appunto, di osservazione. Questo rumore chiaramente compromette il risultato della singola osservazione ma, fortunatamente, è normalmente possibile acquisire molte più osservazioni che incognite ottenendo cosı un sistema sovradimensionato. In queste condizioni, per ottenere una soluzione al problema che minimizzi l'errore, è richiesto l'utilizzo di una tecnica di regressione numerica, per esempio, ai minimi quadrati. In questa prima sezione sono presentate tecniche matematiche ampiamente utilizzate in tutto il libro: per ulteriori dettagli riguardo queste tecniche si può fare riferimento al capitolo 3 incentrato totalmente su questo argomento.

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, detto anche residuo, che questa eventuale soluzione comporterebbe. Non esiste una soluzione generale per un sistema sovradimensionato ma solo soluzioni che minimizzano il residuo sotto una particolare metrica.

Definiamo^1.1 come metrica dell'errore il modulo del residuo

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

La soluzione cosiddetta “ai minimi quadrati” di un sistema lineare (1.1) è rappresenta dal vettore $\mathbf {x}$ che minimizza la distanza euclidea del residuo (1.2) ovvero cercare la soluzione ottima del sistema, nei sensi di una regressione ai minimi quadrati, equivale a trovare il minimo di tale funzione errore al variare di $\mathbf {x}$.

Se si moltiplica l'equazione (1.1) per $\mathbf{A}^{\top}$ si ottiene un sistema lineare “tradizionale” che ammette come soluzione

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

Questo è un primo metodo risolutivo per sistemi sovradimensionati e viene indicato in letteratura come tecnica delle equazioni perpendicolari (normal equations): il problema originale viene ricondotto a un sistema lineare classico dove la matrice dei coefficienti è quadrata e pertanto invertibile con tecniche classiche. Tuttavia la soluzione proposta in equazione (1.3) è numericamente instabile in quanto $\mbox{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.

Si può dimostrare che una soluzione $\mathbf {x}$, meglio condizionata e che minimizza la funzione (1.2), esiste e vale:

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

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.5)

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.

È necessario precisare fin da subito che nel minimizzare la quantità (1.2) non si è fatta nessuna ipotesi sulla distribuzione del rumore all'interno delle varie componenti di cui la matrice è composta: senza tale informazione non c'è garanzia che la soluzione sarà ottima dal punto di vista statistico. Senza ipotesi sulla distribuzione del rumore, la soluzione ottenuta con questa minimizzazione è infatti una soluzione puramente algebrica che minimizza appunto un errore algebrico (algebraic error).

È possibile ottenere una soluzione leggermente migliore dal punto di vista statistico quando il rumore è gaussiano bianco a media nulla e si conosce il valore della varianza del rumore su ogni osservazione. In questo caso è possibile assegnare ad ogni equazione del sistema pesi differenti, moltiplicando ogni riga del sistema per un opportuno peso in modo da pesare in maniera differente ogni dato acquisito. Discussione più approfondita su questo argomento si trova in sezione 3.2 e in generale nel capitolo 2 si affronterà il caso generale dove si conosce il modo con cui l'errore sui dati osservati incide sulla stima dei parametri.

Esistono invece delle tecniche stabili basate su fattorizzazioni che permettono di ricavare la soluzione partendo direttamente dalla matrice $\mathbf{A}$.

Usando per esempio la fattorizzazione QR, algoritmo notoriamente stabile dal punto di vista numerico, della matrice $\mathbf{A}$ il problema originale (1.1) si trasforma nel problema $\mathbf{Q}\mathbf{R} \mathbf{x} = \mathbf{y}$ e la soluzione si può ricavare da $\mathbf{R}\mathbf{x} = \mathbf{Q}^{\top} \mathbf{y}$, sfruttando l'ortogonalità della matrice $\mathbf{Q}$. Nella fattorizzazione QR vige la relazione $\mathbf{R}^{\top}\mathbf{R}=\mathbf{A}^{\top}\mathbf{A}$ ovvero $\mathbf{R}$ è fattorizzazione di Cholesky di $\mathbf{A}^{\top}\mathbf{A}$: attraverso questa relazione si può ricavare infine la pseudoinversa in maniera esplicita.

Attraverso invece la Decomposizione ai Valori Singolari Singular Value Decomposition (SVD), la matrice sovradimensionata $\mathbf{A}$ viene scomposta in 3 matrici dalle proprietà 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 $\mathbf{A} \mathbf{A}^\top$, 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}^{*}$$ (1.6)

dove la pseudoinversa di una matrice diagonale $\mathbf{S}^{+}$ equivale alla sua inversa ovvero una matrice diagonale costituita dai reciproci dei rispettivi valori.

Riassumendo, i modi per risolvere un sistema lineare sovradimensionato sono

• Usando le normal-equations $(A^\top A) x = A^{\top} b$ ovvero la pseudo-inversa di Moore-Penrose $x = (A^\top A)^{-1} A^{\top} b$ (siccome $(A^{\top} A)$ è semidefinita positiva si può usare Cholesky come risolutore);
• La decomposizione QR $A=QR \rightarrow x=R^{-1} Q^{\top} b$;
• la decomposizione SVD $A= U \Sigma V^{\top} \rightarrow x=A^{+} b = V \Sigma^{+} U^{\top} b$;
• La decomposizione di Cholesky $A^{\top} A x = U^{\top} U x = A^{\top} b \rightarrow x = U^{-1}\left( U^{-\top} A^{\top} b \right)$ (essendo $U$ la matrice triangolare superiore $U^{-1}$ e $U^{-\top}$ sono facili da calcolare);

Dal punto di vista del calcolo numerico si può ridurre il numero di condizione della matrice da invertire scalando le colonne di A.

Dettagli ulteriori sulla pseudoinversa di Moore-Penrose possono essere trovati in molti libri, per esempio in (CM09) o nel testo fondamentale di calcolo numerico (GVL96).

Esaminiamo ora il caso in cui il sistema lineare da risolvere sia invece omogeneo.

Un sistema lineare omogeneo ha la forma

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

e normalmente la soluzione ovvia $\mathbf{x}=0$, che è si ottiene anche attraverso l'equazione (1.4), 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, rappresentate un sottospazio vettoriale ovvero il kernel di $\mathbf{A}$. Il vettore generatore del sottospazio è conosciuto a meno di uno o più fattori moltiplicativi (a seconda della dimensione del sottospazio nullo). Per ottenere una soluzione unica è necessario imporre un vincolo aggiuntivo, per esempio $\vert \mathbf{x} \vert = 1$, tale da poter cosı formalizzare

$$\hat{\mathbf{x}} = \argmin_{\mathbf{x}} \Vert \mathbf{A}\mathbf{x} \Vert^{2}$$ (1.8)

con il vincolo $\vert \mathbf{x} \vert = 1$ ovvero realizzando una minimizzazione vincolata.

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

Gli autovettori associati a valori singolari nulli della matrice $\mathbf{S}$ rappresentano pertanto il kernel della matrice stessa e il numero di autovalori nulli rappresenta la dimensione del kernel stesso. Va notato come nell'equazione (1.6) la presenza di zeri nella matrice diagonale $\mathbf{S}$ fosse problematica: ora si capisce che tale presenza è sintomo del fatto che una delle componenti del problema è totalmente incorrelata con la soluzione e, in quanto tale, potrebbe essere trascurata: tale risultato infatti sarà utilizzato nella sezione 2.10.1 nella trattazione dell'agoritmo PCA.

La soluzione del sottospazio di $\mathbf{A}$ è pertanto

$$\mathbf{x} = \sum_{i=1}^{N}{\beta_i \mathbf{v}_i}$$ (1.9)

dove $\mathbf{v}_i$ sono le colonne della matrice $\mathbf{V}$, vettori singolari “destri”, di $\mathbf{A}$ corrispondenti ad $N$ valori singolari, autovalori di 0 della matrice $\mathbf{A}^{\top}\mathbf{A}$.

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

Footnotes

...Definiamo^1.1

la motivazione di questa scelta verrà discussa nei dettagli nel capitolo sullo studio dei modelli.