\chapter{Decomposizione in valori singolari}
 
In questa sezione ci occuperemo di un tipo di fattorizzazione basato sulle proprietà spettrali, e quindi differente
dalle fattorizzazioni viste fino ad ora (come la QR o la LU). Cominceremo esponendo un caso concreto in cui
ci si trova nella necessità di utilizzarla.
 
\section{Problemi lineari ai minimi quadrati}
Consideriamo il seguente problema; sia $A$ una matrice $m \times n$ con $m \geq n$ e $b$ un vettore
di $\R^m$. Vogliamo trovare $x$ tale che $Ax = b$. In generale questo problema è risolubile solo se
$b \in \imm{A}$ vista come applicazione lineare, ma non è conveniente affrontare il problema risolvendo
effettivamente il sistema lineare. Possiamo quindi riformulare la nostra richiesta in questo modo: cerchiamo
$x$ tale che
\[
 x = \min_{x \in \R^n}( || Ax - b ||_2 )
\]
Se $x$ è effettivamente soluzione del sistema lineare avremo $||Ax - b||_2 = 0$, ed in caso contrario
avremo un $x$ tale che $Ax$ sia una buona\footnote{dove buona significa che non si potrebbe trovare
un altro $x'$ tale ceh $Ax'$ sia più vicino (nel senso della norma 2) a $b$ di quanto lo sia già $Ax$}
approssimazione di $b$.
 
Analizziamo alcune tecniche per la risoluzione di questo problema
 
\subsection{Fattorizzazione QR}
Consideriamo la fattorizzazione QR di $A = QR$. Osserviamo che in questo caso $Q \in \matr{\R}{n}{n}$ e
$R \in \matr{\R}{m}{n}$. Osservando che minimizzare la norma 2 è la stessa cosa che minimizzarne il
quadrato si ottiene
\[
 ||Ax - b||_2^2 = ||QRx - b||_2^2 = ||Q(Rx - \trasp Qb)||_2^2 = ||Rx - \trasp Qb||_2^2
\]
dove l'ultimo passaggio è giustificato dal fatto che $Q$ è unitaria.
Consideriamo che la matrice $R$ e il vettore $\trasp Qb$ avranno una struttura di questo tipo:
\[
 R = \left[ \begin{array}{c}
             \hat R \\
             \trasp 0
            \end{array} \right]
\qquad \trasp Qb = \left[ \begin{array}{c}
                            w_1 \\
                            w_2
                          \end{array} \right]
\]
dove $\hat R$ è una matrice quadrata $n \times n$, $w_1$ un vettore con $n$ componenti e $w_2$ un vettore
di $\R^{m-n}$. Possiamo quindi riscrivere la relazione precedente in questo modo
\[
 ||Ax - b||_2^2 = ||\hat Rx - w_1||_2^2 + ||w_2||^2
\]
Abbiamo ricondotto il problema a minimizzare la norma di $\hat Rx - w_1$. Se $\hat R$ è invertibile allora
la soluzione al problema è evidentemente unica ed è la soluzione del sistema lineare $\hat Rx = w_1$.
Osserviamo che $\hat R$ è invertibile se e solo se lo è $\rk{A} = n$; per ora assumeremo che
questo sia vero, riservandoci di analizzare il caso più generale in seguito.
 
In conclusione la fattorizzazione QR ci ha permesso di ricondurre il problema alla risoluzione di un sistema
lineare quadrato di ordine $n$. Analizziamo ora un diverso approccio.
 
\subsection{Sistema delle equazioni normali}
Questo sistema deriva dalla seguente osservazione. Supponiamo di avere $x$ soluzione del sistema $Ax = b$. Allora
$x$ deve essere anche soluzione del sistema quadrato $\trasp AAx = \trasp Ab$, che si ottiene semplicemente moltiplicando
a sinistra la precedente relazione per $\trasp A$.
Osserviamo ancora una volta che scrivendo la fattorizzazione QR di $A$ e supponendo che $A$ abbia rango massimo si
ottiene lo stesso risultato di prima:
\[
 A = QR \qquad \trasp AAx = b \iff \trasp RRx = \trasp R\trasp Qb \iff \trasp R\hat R = \trasp Rw_1 \iff
 \hat Rx = w_1
\]
e quindi si ha l'unicità della soluzione e la risolubilità del sistema lineare. In generale però, il calcolo
di $\trasp AA$ potrebbe esserer oneroso e, soprattutto, nessuno ci garantisce in una situazione generale
che il rango di $A$ sia massimo.
 
\subsection{Teorema di decomposizione in valori singolari}
Possiamo pensare di voler migliorare l'idea delle equazioni normali. In altre parole, vogliamo trovare il
modo di evitare il calcolo di $\trasp AA$. Consideriamo dunque il seguente
\begin{te}[Decomposizione in valori singolari]
Sia $A \in \matr{\R}{m}{n}$. Esistono $U \in \matr{\R}{m}{n}$ e $V \in \mat{\R}{n}$ ortogonali e
$\Sigma \in \matr{\R}{m}{n}$ tali che
\[
 \Sigma = \left[ \begin{array}{cccc}
              \sigma_1 & & & \\
              & \sigma_2 & & \\
              & & \ddots & \\
              & & & \sigma_n \\ \hline
              \multicolumn{4}{c}{\trasp 0} \\
             \end{array} \right]; \qquad \textrm{dove} \: \sigma_i \geq \sigma_j \ \forall i < j \qquad \textrm{e} \qquad
A = U\Sigma\trasp V
\]
\end{te}
 
\begin{os}
 Osservando che $\Sigma$ è quasi diagonale si capisce che in un certo senso questa è una decomposizione spettrale.
 In particolare si può scrivere $\trasp AA$ come
 \[
  \trasp AA = \trasp V\Sigma\trasp U U \Sigma V = \trasp V \left[ \begin{array}{cccc}
              \sigma_1^2 & & & \\
              & \sigma_2^2 & & \\
              & & \ddots & \\
              & & & \sigma_n^2 \\ \hline
              \multicolumn{4}{c}{\trasp 0} \\
             \end{array} \right] V
 \]
e quindi se prendiamo $\lambda_i$ gli autovalori di $\trasp AA$ (che sono reali perché la matrice $\trasp AA$
è simmetrica) ordinati in modo decrescente abbiamo che $\sigma_i = \sqrt{\lambda_i}$ per ogni $i = 1 \ldots n$.
\end{os}
Inoltre questa fattorizzazione mette in luce altre proprietà della matrice $A$. In particolare:
\begin{os}
 Se $A$ ha rango $k$, allora $\sigma_1 \ldots \sigma_k$ sono tutti diversi da $0$ mentre $\sigma_{k+1} \ldots
\sigma_n$ sono $0$. Questo si nota dal fatto che $\trasp AA$ ha lo stesso rango di $A$ ed il suo rango
 è pari al numero\footnote{inteso come numero moltiplicato per la molteplicità algebrica} di autovalori diversi da $0$.
 L'ordinamento che abbiamo richiesto sui $\sigma_i$ ci permette di concludere.
\end{os}
 
Possiamo ora osservare come cambia il problema lineare ai minimi quadrati quando introduciamo questa
nuova fattorizzazione.
\[
 ||Ax - b||_2^2 = ||U\Sigma\trasp Vx - b||_2^2 = ||\Sigma\trasp Vx - \trasp Ub||_2^2
\]
e ponendo $z = \trasp V x$ e $\trasp Ub = \left[ \begin{array}{c} w_1 \\ w_2 \end{array} \right] $ si ottiene
\[
 \Bigg|\Bigg|\Sigma z - \left[ \begin{array}{c} w_1 \\ w_2 \end{array} \right] \Bigg|\Bigg|_2^2 =
 \Bigg|\Bigg|\left[ \begin{array}{ccc} \sigma_1 & & \\ & \ddots & \\ & & \sigma_n \end{array} \right] z - w_1\Bigg|\Bigg|_2^2 + ||w_2||_2^2
\]
Se il rango di $A$ è $n$ allora la matrice $\diag{\sigma_1 \ldots \sigma_n}$ è invertibile e quindi per
minimizzare il primo fattore è sufficiente risolvere un sistema diagonale di ordine $n$.
 
Se però $A$ è di ordine $k < n$ allora si può ripartizionare $\trasp Ub$ in una parte con $k$ componenti
e una parte con le restanti $n-k$. Riscrivendo l'espressione si ottiene che si può minimizzare il primo fattore
risolvendo un sistema di $k$ equazioni e scegliendo arbitrariamente le restanti $n-k$.
\begin{os} Come regola generale si preferisce scegliere queste ultime componenti tutte nulle in modo
 da minimizzare la norma della soluzione. \`E infatti generalmente preferibile ottenere una soluzione di norma piccola.
\end{os}
Dalle relazioni sopra possiamo anche ottenere una espressione esplicita della soluzione $z$. In particolare
\[
 z_j = \frac{w_{1,j}}{\sigma_j} = \frac{\trasp U_j b }{\sigma_j} \quad \text{dove} \ U_j \ \text{è la $j$-esima
colonna di} \ U
\]
Ora possiamo ricordare che $x = Vz$ e quindi se $A$ è di rango massimo
\[
 x = Vz = \sum_{k=1}^{n} \frac{V_kU_kb}{\sigma_k} = \underbrace{\left( \sum_{k=1}^{n} \frac{V_kU_k}{\sigma_k} \right)}_{A^+} b
\]
e chiamiamo $A^+$ \emph{pseudo-inversa} di $A$. Più precisamente
\[
 A^+ = V \Sigma^{-1} \trasp U
\]
Si verifica subito infatti che $AA^+ = I$ e $A^+A = I$\footnote{ma attenzione! Non è vero che $AA^+ = A^+A$ perchè
queste sono addirittura due matrici di ordine diverso. Il loro prodotto dà sempre l'identità, ma non dello stesso
spazio.}.
Questa pseudo-inversa di chiama anche \emph{inversa di Moore-Penrose}.
 
Se $A$ non è di rango massimo possiamo ugualmente definire la pseudo-inversa ``fermando'' la somma
a $k$ dove $k$ è il rango di $A$. Abbiamo dunque
\[
 A^+ = \sum_{j=1}^{k} \frac{V_jU_j}{\sigma_j}
\]
e ancora una volta vale l'eguaglianza $x = A^+b$ dove $x$ è il vettore di rango minimo che risolve (nel modo migliore)
 il problema lineare ai minimi quadrati.
 
\subsection{Interpolazione polinomiale}
Abbiamo visto nel corso di analisi numerica che dati $(x_i,y_i)$ con $i = 0 \ldots n$ è sempre possibile
trovare un polinomio di grado al più $n$ tale che $p(x_i) = y_i$ per tutti gli $i$.
 
Nelle applicazioni concrete si ha però a che fare con grandi moli di dati sperimentali e raramente
il grado del polinomio che si vorrebbe usare per creare un modello matematico ha lo stesso ordine
di grandezza. Ci chiediamo: ``\`E allora possibile trovare il polinomio di grado $m \ll n$ che meglio
approssimi la distribuzione dei nostri dati sperimentali?''; la risposta è sì e si tratta di risolvere
un problema lineare ai minimi quadrati con una matrice $m \times n$.
Indagheremo in seguito questa tecnica per il calcolo del miglior polinomio approssimanete un insieme di dati
sperimentali, perché dovremo superare anche il problema del condizionamente della matrice di Van der Monde
che individua il problema\footnote{Questo era un problema anche nel caso standard con matrice quadrata, e rimarrà
complesso anche nel caso della matrice rettangolare}.
 
\subsection{Principali proprietà della SVD} \label{subsec:svdprop}
Abbiamo introdotto la SVD parlando di sistemi lineari perché questo è il contesto in cui è stata originariamente
concepita. Possiamo però osservare che essa ha anche delle interessanti proprietà in molti altri ambiti.
 
Consideriamo la matrice $A \in \matr{\R}{m}{n}$ che individua univocamente un'applicazione lineare
\[
T_A : \R^n \longto \R^m
\]
La \svd\ ci permette di conoscere meglio questa applicazione. Sappiamo che se $A$ è di rango $k$ allora
$A = U\Sigma\trasp V$ e i $\sigma_i$ della matrice $\Sigma$ sono nulli a partire dal ($k+1$)-esimo.
\`E chiaro dunque che la scomposizione di permette di determinare il rango di $A$ e quindi di capire
che l'immagine di $T_A$ avrà dimensione $k$ e il kernel dell'applicazione $n-k$. Ricordando però che
\[
 Ax = U\Sigma\trasp Vx = \sum_{i=1}^{k} \sigma_i U_i\trasp V_i x
\]
possiamo estrarre altre informazioni. Dato che per ogni $x$ $V_i x$ è uno scalare l'immagine di $T_A$
deve essere contenuta nello $\Span{U_1 \ldots U_k}$. Dato che quest'ultimo ha dimensione $k$, che è la
stessa di $\imm{T_A}$, devono coincidere e quindi le prime $k$ colonne di $U$ sono una base per l'immagine.
Analogamente dall'ortogonalità delle colonne di $V$ si ha che se $x = V_j$ con $j > k$ allora $Ax = 0$
e quindi lo $\Span{V_{k+1} \ldots V_{n}}$ è contenuto nel kernel di $T_A$. Ancora una volta per motivi
dimensionali devono coincidere.
 
In sintesi si ha che se $A = U\Sigma\trasp V$ è una matrice di rango $k$ allora
\begin{itemize}
 \item $\Ker{A} = \Span{V_{k+1} \ldots V_n}$;
 \item $\imm{A} = \Span{U_1 \ldots U_k}$;
\end{itemize}
 
Introduciamo ora un'estensione del concetto di norma indotta visto nel corso di Analisi Numerica
\begin{de}
 Sia $A$ una matrice $m \times n$. Si dice \emph{norma di $A$ indotta dalla norma vettoriale $||\cdot||$}
 la seguente funzione
 \[
  \begin{array}{rcl}
||\cdot|| : \matr{\R}{m}{n} &\longto& \R \\
 A & \longmapsto & \displaystyle  \max_{||x|| = 1} ||Ax||
              \end{array}
 \]
\end{de}
Questa è effettivamente una norma, ma la verifica viene qui omessa.
%TODO: Magari si potrebbe anche fare, suvvia!
 
\begin{os}
 Osserviamo la seguente catena di uguaglianze
 \[
 ||Ax||_2 = ||U\Sigma \trasp V x||_2 = ||\Sigma \trasp V x||_2 = ||\Sigma z||_2
\]
Queste ci dicono che, data l'invertibilità e l'ortogonalità di $V$, $||A||_2 = ||\Sigma||_2 = \sigma_1$.
\end{os}
Date queste considerazioni possiamo introdurre anche una generalizzazione del condizionamento
ridefinendole con la nuova norma e con la pseudo inversa. Nel caso della norma due si ottiene che
\[
 \cond{A} = ||A^+|| \cdot ||A||_2 = \frac{\sigma_1}{\sigma_k}
\]
dove $\sigma_k$ è il più piccolo $\sigma_i$ non nullo, ovvero, come al solito, $k$ è il rango di $A$.
Un'analisi dell'errore effettuata sul problema lineare ai minimi quadrati mostrerebbe come effettivamente
l'errore generato dal procedimento dipende proprio da questo coefficiente.
 
\subsection{\tsvd, pulizia del rumore e compressione di immagini}
Le considerazioni fatte nella Sezione~\ref{subsec:svdprop} ci permettono di introdurre delle tecniche
per ridurre il condizionamento di alcuni problemi.
 
Supponiamo di avere ad esempio un sistema lineare $Ax = b$ con $A$ una matrice con un condizionamento
molto alto. Siamo consci di non avere speranza di risolvere in maniera esatta il sistema senza aumentare
la precisione dei tipi che utilizziamo negli algoritmi.
Se ad esempio il condizionamento è nell'ordine di $10^{12}$ e la precisione di macchina di $10^{16}$ sappiamo
che probabilmente avremo circa $4$ cifre decimali corrette (nella migliore delle ipotesi). La soluzione di aumentare
la precisione di macchina (ad esempio portarla a $10^{-32}$) potrebbe essere una soluzione ma appesantirebbe
molto l'esecuzione dell'algoritmo.
 
Consideriamo dunque la \svd di $A = U\Sigma\trasp V$.
Possiamo supporre che $\sigma_n$ sia molto piccolo e che
questo causi il cattivo condizionamento (ricordiamo che $\cond{A} = \frac{\sigma_1}{\sigma_n}$). Possiamo
pensare di troncare la matrice $\Sigma$ sostituendola con una matrice $\hat \Sigma$ identica ad eslusione del posto
$(n,n)$ dove rimpiazziamo $\sigma_n$ con $0$, e risolvere il problema $\hat Ax = U\hat \Sigma\trasp Vx = b$.
Questo \textbf{è un altro problema}, ma potrebbe essere molto meglio condizionato del precedente (ad esempio
se $\sigma_{n-1} \gg \sigma_n$) e talvolta è la scelta migliore.
 
Questa pratica di ``troncamento'' viene chiamata \tsvd, ovvero \svd\ con troncamento. \`E tipicamente
utile nella ricostruzione di immagini (che vedremo più avanti) perché non ci interessa la soluzione esatta
del problema, ma la soluzione esatta di un problema ``semplice'' che si avvicini ai nostri dati.
Questo procedimento viene anche detto di soglia, o di filtro, perché è possbile decidere una soglia minima
sotto la quale i $\sigma_i$ non possono scendere, pena essere messi a $0$. In questo modo, fissata la
soglia $s$, avremo che il condizionamento sarà maggiorato da $\sigma_1s^{-1}$.
 
In generale non è un problema banale scegliere il livello di soglia $s$. Nei casi più semplici si ha
un salto dei valori singolari, che passano da un ordine di grandezza ad uno (o a più) inferiori, dividendo
in due zone ben distinte quelli prossimi a $0$ e quelli $\gg 0$.
 
Osservaiamo che in generale non abbiamo alcuna garanzia che i risultati ottenuti tramite
questa tecnica abbiano un riscontro con il problema che aveva generato i dati di partenza, e va quindi
usata con molta attenzione.
 
Un altro esempio pratico di utilizzo della \tsvd\ è quello della compressione delle immagini. Supponiamo di avere
un'immagine rappresentata come una matrice di pixel $A$\footnote{non esiste un motivo per cui la matrice dovrebbe
essere quadrata, è solo per rendere più chiaro l'esempio} di dimensione $n \times n$. Possiamo scrivere
$A = U\Sigma\trasp V$ e pensare di limitarci ad una matrice di rango $k$ e quindi tronvare tutti i $\sigma_i$
con $i > k$. In questo modo potremo rappresentare l'immagine considerando solo le prime $k$ colonne di $U$, le
prime $k$ righe di $V$ e i $\sigma_1 \ldots \sigma_k$.
Per fare un esempio pratico supponiamo $n = 1024$, ovvero un'immagine da 1 Megapixel. Questa occuperebbe in memoria
(suppoendo che sia in bianco e nero a 256 colori) $1 MB$. Decidendo di comprimerla con una matrice di rango $15$
avremmo invece una dimensione di $15 KB$! Ovviamente l'immagine risultante darebbe solamente un'idea di quella originale,
ma si potrebbe pensare di trovare un punto d'incontro fra dimensione e qualità\footnote{e magari di raffinare anche
il metodo di compressione, ma questo lo vedremo in seguito}.