\chapter{Calcolo degli autovalori di matrici strutturate}
 
Abbiamo visto nell'analisi dei sistemi lineari vari metodi per utilizzare la struttura
delle matrici per ottimizzare il procedimento e diminuire il costo computazionale. \\
Questo in generale è più difficile nel calcolo degli autovalori. L'unica struttura che siamo finora
riusciti a sfruttare per abbassare il costo è stata la tridiagonalità.
 
\section{Strutture note}
In generale si parla di matrice strutturata quando una matrice è individuata da un numero di
parametri dell'ordine minore di $n^2$. \\
Ricapitoliamo le strutture principali viste fino ad ora e i relativi vantaggi nella risoluzione
di sistemi lineari.
\begin{description}
 \item[Matrici sparse] Parliamo di matrici sparse quando ci sono pochi elementi non nulli, ovvero un numero
 con un ordine di grandezza strettamente inferiore a $n^2$. In questo caso riusciamo spesso ad ottimizzare
 i metodi di moltiplicazione matrice per vettore e quindi possiamo trarre vantaggio dall'applicazione
 di metodi iterativi;
 
 \item[Van der Monde] Non abbiamo analizzato in dettaglio\footnote{ad eccezione delle matrici di Fourier per il calcolo
 della \dft, dove si riesce a scendere ad un costo $O(n\log n)$.} queste matrici ma esiste un metodo (tramite un'opportuna
 rappresentazione dell'inversa) per risolvere un sistema lineare con un costo $O(n^2)$.
 
 \item[Toeplitz o H\"ankel] Abbiamo visto queste matrici nella Sezione~\ref{subsec:toeplitz} e tramite la
 trasformata discreta di Fourier abbiamo mostrato che il costo del prodotto matrice vettore è molto basso.
 Possiamo quindi ipotizzare di applicare dei metodi iterativi per risolvere il sistema.
 
 \item[Diagonali con correzione di rango 1] Queste martici sono piuttosto interessanti e sono le prime di
 cui ci occuperemo. Sono appunto composte dalla somma di una matrice diagonale $D$ con una matrice di rango $1$
 $u\trasp u$;
\end{description}
 
\section{Il metodo di Lanczos}
\subsection{Un'alternativa al metodo di Householder}
Introduciamo ora un metodo alternativo a quello di Householder per la tridiagonalizzazione di una
matrice simmetrica. \\
Osserviamo che se $A$ è una matrice simmetrica in $\mat{\R}{n}$, allora esiste una matrice $T \in \mat{\R}{n}$
tridiagonale simmetrica e una matrice unitaria $Q$ tali che
\[
 A = QT\trasp Q
\]
e quindi anche
\[
 AQ = QT
\]
Osserviamo ora solo la prima colonna di quest'uguaglianza, ovvero $AQe_1 = Aq_1 = QTe_1$. Se scriviamo
\[
 T = \left[ \begin{array}{cccc}
             \alpha_1 & \beta_1 & & \\
             \beta_1  & \ddots  & \ddots & \\
              & \ddots & \ddots & \beta_{n-1} \\
              & & \beta_{n-1} & \alpha_n \\
            \end{array} \right]
\]
Possiamo riscrivere la relazione come $Aq_1 = \alpha_1 q_1 + \beta_1 q_2$ e, una volta scelto $q_1$ vettore
di norma $1$,
si ha che $\trasp q_1 A q_1 = \alpha_1 \trasp q_1 q_1 + \beta_1 \trasp q_1 q_2 = \alpha_1$; possiamo quindi
calcolare $\alpha_1$ da cui possiamo poi ricavare $\beta_1$ considerando che $\beta q_2 = Aq_1 - \alpha_1 q_1$
e quindi
\[
\beta_1 = ||Aq_1 - \alpha_1 q_1||_2 \qquad \text{e} \qquad q_2 = \frac{Aq_1 - \alpha_1 q_1}{\beta_1}
\]
Ripetendo il procedimento con la seconda colonna (ora che conosciamo $q_2$) e poi continuando ancora si ottiene
la regola generale
\[
 \alpha_i = \trasp q_i A q_i \qquad \beta_i = ||Aq_i - \alpha_i q_i||_2 \qquad q_{i+1} = \frac{Aq_i - \alpha_i q_i}{\beta_i}
\]
e si può quindi ricostruire tutta la matrice $Q$ e la matrice $T$.
 
Il costo computazionale dominante nello svolgimento di queste operazioni è dato dai prodotti matrice vettore
che costano generalmente $O(n^2)$ e portano quindi ad un costo complessivo del metodo di $O(n^3)$ (lo stesso
costo del metodo di Householder). \\
Questo metodo non viene però generalemente utilizzato per calcolare la tridiagonalizzazione. Il motivo è che
non è numericamente stabile, ed in generale la matrice $Q$ ottenuta non è ortogonale.
 
\begin{os}
 A questo punto è naturale chiedersi che rapporto esiste fra la matrice calcolata con il metodo di Householder
 e quella calcolata con il metodo di Lanczos, al variare del primo vettore $q_1$ scelto all'inizio. La risposta
 che è possibile verificare è che scegliendo $q_1 = e_1$ le matrici differiscono per una matrice di fase.
\end{os}
 
Nonostante in generale il metodo non venga usato per la sua poca precisione esiste un particolare caso in cui
risulta essere utile, ed è precisamente il seguente
 
\subsection{Il quoziente di Rayleigh}
Cominciamo con una definizione
\begin{de}
 Si dice \emph{quoziente di Rayleigh di $A$ e $x$} e si indica con $\ray{A}{x}$ lo scalare
 \[
  \ray{A}{x} = \frac{\trasp x A x}{\trasp xx}
 \]
\end{de}
Si può osservare che il minimo del quoziente di Rayleigh su tutto $\R^n$ corrisponde al modulo dell'autovalore minimo
e il massimo al modulo dell'autovalore massimo. \\
Osserviamo poi che preso un generico sottospazio $k$-dimensionale $S \subseteq \R^n$ abbiamo che se
\[
 \lambda = \min_{x \in S} \ray{A}{x}
\]
allora $\lambda$ è l'autovalore di modulo minimo su un sottospazio di $\R^n$. Se $S$ è sufficientemente grande
possiamo pensare di usare $\lambda$ come approssimazione di $\lambda_{min}$. \\
Consideriamo $\{ q_1 , \ldots, q_k \}$ base del sottospazio $S$ e la matrice $Q$ con i vettori $q_i$ come colonne.
Se prendiamo $x \in S$ allora
\[
x = \sum_{i=1}^{k} \alpha_i q_i
\]
e quindi se $\alpha = (\alpha_1, \ldots, \alpha_k)$ si ha $x = Q\alpha$.
\[
\frac{\trasp x A x}{\trasp xx} = \frac{\trasp \alpha \trasp Q A Q \alpha}{\trasp \alpha \trasp Q Q \alpha} =
 \frac{\trasp \alpha \trasp Q A Q \alpha}{\trasp \alpha \alpha}
\]
In conclusione è sufficiente minimizzare $\ray{\trasp Q A Q}{x}$ su $\R^k$, ed essendo $\trasp Q A Q$ una matrice
$k \times k$ questo procedimento può risultare sensibilmente più economico rispetto all'idea originale. \\
\begin{de}
 Sia data una matrice $A$ a coefficienti reali\footnote{definiamo il procedimento sui numeri reali solo per non appesantire
 la notazione, ma non c'è nessuna restrizione ad usarlo sui complessi.} e $x$ un vettore di $\R^n$.
 Si dice \emph{sottospazio di Krylov} di $A$ e $v$ di ordine $j$ e si indica con $\kryl{A}{v}{j}$ il sottospazio
 \[
  S = \Span{v, Av, \ldots, A^{j-1}v}
 \]
\end{de}
Osserviamo che in generale $\dim{S} \leq j$. \\
Siamo ora interessati ad utilizzare un sottospazio di Krylov come sottospazio $S$ con cui approssimare il nostro autovalore.
Per farlo però dobbiamo trovare una base ortonormale di $S$; possiamo procedere calcolando la scomposizione \qr\ della
matrice dei vettori che generano:
\[
 \left[ \begin{array}{c|c|c|c}
         \multirow{4}{*}{$v$} & \multirow{4}{*}{$Av$} & \multirow{4}{*}{$\ldots$} & \multirow{4}{*}{$A^{j-1}v$} \\
  & & & \\
  & & &  \\
  & & & \\
        \end{array} \right] = QR
\]
Si può mostrare che le prime $j$ colonne di $Q$ sono una base dello spazio $S$. Per farlo è sufficiente osservare
che essendo $R$ triangolare superiore e $n \times j$ con $j < n$ deve avere le ultime $n - j$ righe nulle. \\
Calcolando ora $\trasp Q A Q = (b_{ij})$ si ottiene un matrice $k\times k$ simmetrica. Osserviamo in particolare
che se $i < j+1$ allora $b_{ij} = \trasp q_i A q_j$ e $Aq_j \in \Span{v, Av, \ldots , A^{j-1}v}$ per costruzione.
In particolare $Aq_j$ è ortogonale a $q_i$ e quindi $b_{ij} = 0$. La simmetria della matrice ci permette
di concludere che $B$ è tridiagonale, e si può verificare che è la matrice $k\times k$ generata dal metodo
di Lanczos dopo $k-1$ iterazioni. \\
A questo punto si può supporre di avere una buona approssimazione dell'autovalore di modulo massimo (anche se
non c'è nessuna garanzia di averla davvero). Questa è di fatto l'unica applicazione del metodo di Lanczos,
data la sua instabilità numerica.
 
\section{Matrici con struttura di rango}\label{sec:struttura_di_rango}
\subsection{Qualche esperimento con il metodo \qr}
Abbiamo spesso incontrato durante il corso matrici diagonali con correzioni di rango 1,
ovvero della forma $D + u\trasp u$. Queste
sono molto interessanti perché si ritrovano in molti problemi computazionali e perché sono inverse
di matrici tridiagonali. \\
Ci domandiamo ora se esiste un metodo efficiente per applicare il Metodo~\qr\ a queste matrici.
Sappiamo di poterle tridiagonalizzare con matrici di Householder con un costo computazionale $O(n^2)$, ma
vogliamo provare ad applicare il \qr\ alla matrice piena senza dover necessariamente passare
per il procedimento di tridiagonalizzazione\footnote{Abbiamo fatto notare nella Sezione~\ref{subsec:qr_costo} che
in generale non è conveniente applicare il metodo ad una matrice piena; cercheremo però di trovare
qualche via più breve per questa particolare classe di matrici.}. \\
Possiamo osservare sperimentalmente che facendo qualche passo del metodo \qr\ (anche utilizzando lo shift)
la struttura non viene completamente persa.
 
Indichiamo con $\alpha_k$ il rango massimo delle sottomatrici quadrate strettamente contenute nella parte
inferiore della matrice della $k$-esima iterazione del \qr, e con $\beta_k$, analogamente, il rango massimo
di quelle superiori.
 
\textbf{Esperimento 1}: Supponiamo di scegliere una matrice qualsiasi (complessa) e di applicare
il metodo. Osserveremo che $\alpha_0 = \alpha_1 = \ldots = \alpha_k = 1$
mentre $\beta_k$ cresce.
 
\textbf{Esperimento 2}: Supponiamo ora di restringere la scelta ad una matrice con $D$ reale, ovvero
$A = D + u\trasp u$ dove $D = \diag{\gamma_1, \ldots, \gamma_n}$ e i $\gamma_i$ sono tutti reali. Allora
otterremo sperimentalmente $\alpha_i = 1$ e $\beta_i = 3$ per ogni $i \geq 2$.
 
\textbf{Esperimento 3}: Analogamente a prima scegliamo $D$ complessa ma con elementi di modulo $1$\footnote{ovvero
quella che abbiamo chiamato \emph{matrice di fase}.}. Otterremo ancora una volta $\alpha_i = 1$ e $\beta_i = 3 \:
\forall i \geq 2$.
 
Come mai i $\beta_k$ in questi ultimi due casi non crescono più di $3$? Cerchiamo di rispondere analizzando separatamente
i vari casi.
 
\subsection{Conservazione del rango nella parte triangolare inferiore}
Ricordiamo dall'Osservazione~\ref{os:qr_simili_ak} che le matrici $A_k$ generate dal metodo \qr\ sono
simili per trasformazione ortogonale. Se assumiamo che $R_k$ sia invertibile possiamo anche mostrare
qualcosa di più
\[
 A_{k+1} = R_k Q_k + \alpha_k I = R_k Q_k R_k R_k^{-1} + \alpha_k R_k R_k^{-1} = R_k A_k R_k^{-1}
\]
e quindi $A_{k+1}$ è simile ad $A_k$ tramite una trasformazione con una matrice triangolare superiore
$R_k$. Questo ci assicura che nella parte inferiore della matrice $A_{k+1}$ il rango venga conservato
(nel senso degli $\alpha_k = 1$).
 
%% TODO: Spiegare perché il rango si conserva, da fare quando anche io lo avrò capito.
 
\subsection{Limitazione del rango nella parte triangolare superiore}
Ci occuperemo ora di spiegare perché il rango di tutte le sottomatrici strettamente contenute nella
parte triangolare superiore è limitato superiormente da $3$, sia nel caso in cui $D$ sia una matrice
reale, sia in quello in cui sia una matrice di fase. Consideriamo il primo e osserviamo cosa succede
al primo passo del metodo. Poniamo $A_0 = D + u\trasp v$; si ha che:
\[
 A_1 = \herm{Q_1} A_0 Q_1 = \herm{Q_1}(D+ u\trasp v)Q_1 = \herm{Q_1} DQ_1 + u_1\trasp{v_1}
\]
Ricordando che $D$ è reale si può concludere che $\herm{Q_1}DQ_1$ è hermitiana e diagonale;
possiamo quindi scrivere il $k$-esimo passo del metodo in questo modo
\[
 A_{k+1} = H_{k+1} + u_{k+1}\trasp v_{k+1}
\]
dove $H_{k+1}$ è una matrice hermitiana. In particolare si ha $H_{k+1} = A_{k+1} - u_{k+1}\trasp v_{k+1}$.
Da quanto visto prima sappiamo che il rango nella parte inferiore di $H_{k}$ non può superare $2$.
Ricordando che $H_k$ è hermitiana si può concludere che il rango è al massimo $2$ nella parte superiore.
Osservando nuovamente che $A_k = H_k + u_k\trasp v_k$ si ottiene che il rango di $A_k$ nella parte superiore
non può superare $3$ e si ha quindi la tesi.
 
La stessa cosa si può mostrare che quando la matrice diagonale $D$ è una matrice di fase. Ci sono due
procedure possibile per effettuare la dimostrazione; una è scrivere $D = QR$ ed osservare che data
l'hermitianità di $D$ si ottiene che $R$ deve forzatamente essere una matrice di fase. Questa via
richiede però tediosi passaggi molto contosi che quindi non svolgeremo. Un'alternativa più elegante
è utilizzare il \emph{Nullity theorem}, dimostrato attorno al 1960. Questo però non fa parte del
programma del corso e quindi non seguiremo neanche questa via.
 
\subsection{Metodo \qr\ per il calcolo delle radici dei polinomi}
Il motivo principale (e storicamente il primo) per cui ci si è interessati a queste matrici
è stato la ricerca delle radici dei polinomi tramite l'uso delle matrici Companion. \\
Le matrici Companion hanno una struttura come la seguente:
\[
F = \left[ \begin{array}{cccc}
            0 & \cdots & \cdots & \times \\
            1 & \ddots & & \vdots \\
	       &\ddots & 0 & \times \\
              & & 1 & \times
           \end{array} \right]
\]
La matrice $F$ si può scrivere come una matrice unitaria più una correzione di rango $1$ nel seguente modo
\[
F = Q + u\trasp e_n = \left[ \begin{array}{cccc}
                              0 & & & 1 \\
                              1 & \ddots & & \\
                               & \ddots & \ddots & \\
                               & & 1 & 0 \\
                             \end{array} \right]
+ \left[ \begin{array}{cccc}
                               \: &\:  &\:  & v_1 \\
                               &  & & \vdots \\
                               &  &  & \vdots \\
                               & &  &  v_n \\
                             \end{array} \right]
\]
e quindi sarebbe conveniente poter applicare delle ottimizzazioni al \qr\ sfruttando le osservazioni fatte
sul rango delle matrici contenute nella parte inferiore e superiore della matrice che viene iterata.
 
Supponendo di voler determinare le radici del polinomio $p(x) = p_0 + p_1x + \ldots + p_nx^n$ si ottiene
la matrice
\[
F = \left[ \begin{array}{cccc}
            0 & \cdots & \cdots & \frac{-p_0}{p_n}\\
            1 & \ddots & & \vdots \\
	       &\ddots & 0 & \vdots \\
              & & 1 & \frac{p_{n-1}}{p_n}
           \end{array} \right]
\]
In generale, affrontando il problema della ricerca degli autovalori, usiamo il teorema di Bauer-Fike (Teorema~\ref{te:BauerFike})
per assicurarci di ottenere il risultato con una buona approssimazione. Questo però non è sempre sufficiente
se $p_n$ è piccolo, perché la norma della matrice può crescere arbitrariamente. Come affrontare questo problema?
 
Una soluzione può essere evitare di ricondursi ad un problema agli autovalori analizzando il seguente \emph{problema
generalizzato agli autovalori}\footnote{In generale si dice \emph{problema generalizzato agli autovalori} un problema
del tipo $\deter{A - \lambda B} = 0$. Se $B$ è invertibile questo può sempre essere ricondotto a trovare gli
autovalori di $B^{-1}A$, ma non sempre è conveniente. }:
\[
 \deter{F - \lambda I} = 0 \iff \deter{p_n F - \lambda( I + p_n e_n\trasp e_n)} = 0
\]
Parleremo dei possibili metodi per analizzare questo problema nel Capitolo~\ref{cap:autovalori_generalizzato}.