\chapter{Decomposizione in valori singolari} \label{ch:svd}
 
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 essere 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]
Siano $m \geq n$ e $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 \\
             \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 \\
             \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_k\trasp U_kb}{\sigma_k} = \underbrace{\left( \sum_{k=1}^{n} \frac{V_k\trasp U_k}{\sigma_k} \right)}_{A^+} b
\]
e chiamiamo $A^+$ \emph{pseudo-inversa} di $A$. Più precisamente
\[
 A^+ = V \Sigma^{-1} \trasp U
\]
% TODO: Controllare che le uguaglianze scritte sull'inversa abbiano senso. Una è giusta di sicuro,
% ma quale?
Si verifica subito infatti che 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_j\trasp U_j}{\sigma_j}
\]
e ancora una volta vale l'eguaglianza $x = A^+b$ dove $x$ è il vettore di rango minimo che risolve (nel senso
che minimizza il valore della norma)
 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
(supponendo 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} (si veda ad esempio la Figura~\ref{fig:camosci}).
 
\begin{figure}[ht!]
\begin{center}
  \includegraphics[scale=0.2]{camosci/camoscio_100.png}
  \includegraphics[scale=0.2]{camosci/camoscio_50.png}
  \includegraphics[scale=0.2]{camosci/camoscio_25.png}
  \includegraphics[scale=0.2]{camosci/camoscio_15.png}
  \includegraphics[scale=0.2]{camosci/camoscio_8.png}
 \caption{Una successione di camosci di 300 pixel di lato e di rango rispettivamente $100$, $50$, $25$, $15$ e $8$; le
 immagini sono state ridimensionate per stare nella pagina, e sono state realizzate con un programma
 identico a quello realizzato in laboratorio durante il corso.} \label{fig:camosci}
\end{center}
\end{figure}
 
 
\section{Calcolo della decomposizione in valori singolari}
In questa sezione introdurremo dei metodi per il calcolo della \svd\ e faremo qualche altra considerazione
sui suoi possibili utilizzi.
 
\subsection{Approssimazione di rango fissato}
Supponiamo di avere una matrice $A \in \mat{\R}{n}$ e di volerne trovare una approssimazione efficiente.
Un'idea può essere quella di rappresentarla con una matrice di rango fissato $k < n$ che le assomigli il
più possibile. \\
Consideriamo la procedura seguita anche alla fine della sezione precedente
\begin{enumerate}
 \item Calcoliamo la fattorizzazione \svd\ di $A = U \Sigma V$;
 \item Tronchiamo la $\Sigma$ ponendo $\sigma_{k+1} \ldots \sigma{n} = 0$ e chiamiamo $\tilde \Sigma$
 la nuova matrice;
 \item Approssimiamo $A$ con $U \tilde \Sigma V = \sum_{i=1}^{k} \sigma_j U_j \trasp{V_j}$;
\end{enumerate}
\begin{os}
 Questa approssimazione di $A$ è molto pratica nel caso non si abbia sufficiente spazio per
 salvare $A$ nella memoria di un calcolatore. Suggerisce infatti un metodo immediato per ridurre
 (circa) di un fattore $\frac{n}{k}$ lo spazio occupato. \`E infatti sufficiente salvare le prime $k$
 colonne di $U$, le prime $k$ righe di $A$ e i primi $k$ valori sulla diagonale di $\Sigma$.
\end{os}
Una domanda naturale a questo punto è: ``Quanto differiscono le due matrici?''. Osserviamo che
\[
 ||A - B||_2 = ||U\Sigma \trasp V - U \tilde \Sigma \trasp V||_2 = ||U[\Sigma - \tilde \Sigma] \trasp V||_2
 = \sigma_{k+1}
\]
Se $\sigma_{k+1}$ è sufficientemente piccolo abbiamo quindi un piccolo errore. Questo, di fatto, si
verifica in molte applicazioni pratiche dove i $\sigma_i$ piccoli possono rappresentare i disturbi
all'informazione che si intendeva analizzare (come nel caso della trasmissione di segnali).
 
In realtà possiamo dire di più. Sia $B$ una matrice di rango $k$. Sappiamo che
\[
 \Dim{\Ker{B}} + \Dim{\imm{B}} = n \qquad \Dim{\Ker{B}} \geq n - k
\]
Considero il sottospazio $S = \Span{V_1,\ldots,V_{k+1}} \subseteq \R^n$. Per una questione di dimensioni
deve esistere un vettore non nullo $z \in S \cap \Ker{B}$, e possiamo supporre $||z||_2 = 1$.
Osserviamo che
\[
|| (A - B)z ||_2 = ||Az||_2 = ||\sum_{i=0}^{k+1} \sigma_i U_i \trasp V_i z||_2 \geq
 || \sigma_{k+1} U_i \trasp V_i z ||_2 = \sigma_{k+1}
\]
In particolare abbiamo mostrato che la nostra approssimazione era la migliore possibile, nel senso che
per ogni altra approssimazione ha una distanza maggiore o uguale da $A$ (secondo la norma 2).
 
\subsection{Metodi di calcolo, errori e costo computazionale}
Abbiamo visto che la \svd\ è una sorta di decomposizione spettrale della matrice $A$. Più precisamente
i $\sigma_i$ che si trovano sulla diagonale di $\Sigma$ sono le radici degli autovalori di $\trasp AA$
(che è semidefinita positiva). \\
Un primo modo per affrontare il calcolo della fattorizzazione potrebbe quindi essere
\begin{enumerate}[(a)]
 \item Calcolare $\trasp AA$;
 \item Determinare la decomposizione spettrale di $\trasp AA = V D \trasp V$;
 \item Calcolare la fattorizzazione $QR$ di $AV = UR$;
\end{enumerate}
Osserviamo ora che $\trasp R\trasp U U R = \trasp RR = \trasp V \trasp AA V = D$; $\trasp RR$ è quindi
la fattorizzazione di Cholesky di $D$. Si può mostrare che la diagonalità di $D$ implica che anche $R$ sia
diagonale e quindi $A = UR\trasp V$ è quasi la \svd\ di $A$. In effetti non abbiamo ancora considerato
il problema dell'ordine (vorremmo $\sigma_i \geq \sigma_j$ se $i \geq j$).
Questo problema sarebbe risolto nel momento in cui la colonna di $AV$ di norma maggiore fosse in prima posizione
(e conseguentemente le altre in ordine di norma discendente).
Possiamo quindi introdurre una opportuna matrice di permutazione $P$ e riscrivere la relazione
nel seguente modo
\[
 AVP = U\Sigma
\]
In questo modo la prima colonna di $AVP$ è quella di norma maggiore e quindi i $\sigma_i$ saranno
in ordine decrescente.
 
\begin{os}
 Per eseguire il procedimento sopra descritto abbiamo bisogno di conoscere la matrice $\trasp AA$, che tipicamente
 può essere grande (ad esempio se $m \gg n$). Calcolare la sua decomposizione spettrale potrebbe introdurre
 molti problemi di stabilità, e occupare molto spazio in memoria. Sarebbe quindi interessante trovare un
 metodo per avere una decomposizione spettrale di $\trasp AA$ senza doverla effettivamente calcolare.
\end{os}
 
Supponiamo ora di avere $U,V$ matrici ortogonali e $B$ matrice bidiagonale tali che
\[
 A = UB\trasp V
\]
allora si ha che
\[
 \trasp AA = \trasp V \trasp B B V
\]
e $\trasp B B$ è una matrice tridiagonale; avere la decomposizione spettrale di $\trasp B B$ è quindi equivalente
ad avere la decomposizione spettrale di $\trasp A A$!
 
Sappiamo inoltre che calcolare la decomposizione spettrale di una matrice tridiagonale simmetrica è particolarmente
semplice.
 
L'unico pezzo che manca è mostrare come è possibile calcolare le matrici $U$ e $V$ che bidiagonalizzano $A$.
Possiamo usare le matrici di Householder in questo modo:
\begin{enumerate}[(a)]
 \item Calcoliamo $P_1$ tale che $P_1 Ae_1 = \alpha e_1$, e quindi $P_1A$ abbia la prima colonna
 nulla ad eccezione dell'elemento in posizione $(1,1)$;
 \item Calcoliamo $Q_1$ tale che $P_1 A Q_1$ abbia la prima riga nulla ad esclusione dei primi
 due elementi (in posizione $(1,1)$ e $(1,2)$);
 \item Iteriamo il procedimento sopra lavorando sulla seconda colonna e sulla seconda riga, e poi di seguito
 per tutte le righe e colonne;
\end{enumerate}
Si verifica che il procedimento sopra porta effettivamente ad una matrice bidiagonale e che $\trasp BB$
è tridiagonale simmetrica.
 
Analizziamo ora il costo computazionale di questo procedimento. Il primo passo consiste nel calcolo delle
due matrici ortogonali che bidiagonalizzano $A$. Il costo di ogni passo del procedimento è il calcolo
di due matrici di Householder e quindi $O(n)$. Viene ripetuto $n$ volte fino ado ottenere un
costo complessivo di $O(n^2)$. Diventa poi $O(mn^2)$ perchè bisogna effettuare le moltiplicazioni.
 
Analogamente il costo del calcolo degli autovalori con il metodo \qr\ è $O(n^2)$; il metodo richiederebbe
però di calcolare la matrice $\trasp BB$ che potrebbe essere piuttosto grande. In realtà esiste un metodo
(che qui non esponiamo) per applicare il metodo \qr\ a $\trasp BB$ senza calcolarla esplicitamente ma
conoscendo solamente $B$, lasciando il costo totale a $O(mn^2)$.
 
Un problema da affrontare, infine, è quello dell'errore. Il teorema di Bauer-Fike (Teorema~\ref{te:BauerFike}) ci
assicura infatti che il problema del calcolo degli autovalori di una matrice simmetrica è ben condizionato,
ma assicurandoci solo una maggiorazione dell'errore assoluto, e non una dell'errore relativo. Sfortunatamente
nelle applicazioni si è spesso interessati a calcolare con precisione i valori singolari piccoli e questo
potrebbe essere un problema. Negli ultimi anni sono state sviluppate tecniche per il calcolo della \svd\ che
hanno tenuto conto di questo problema ma non sono ancora state implementate in programmi come \matlab\ o librerie
come \lapack.