\chapter{Il calcolo degli autovalori}
 
In questo capitolo ci occuperemo di presentare i metodi principalmente usati nel calcolo
effettivo degli autovalori di matrici Hermitiane tridiagonali e in forma di Hessenberg.
 
\section{Il metodo di Sturm}
Introdurremo ora un metodo efficiente per il calcolo di singoli autovalori di una matrice
hermitiana in forma tridiagonale. Per tutta la sezione tratteremo il caso di una matrice
reale per semplicità, ma il caso complesso è risolubile in modo analogo.
 
Abbiamo visto nella Sezione~\ref{sec:eqsecolare} un metodo efficiente per valutare il polinomio
caratteristico di una matrice di questo tipo con costo $O(n)$. Possiamo affinare la nostra osservazione
notando che anche la valutazione di tutti i polinomi $p_k$ in un punto $x$ fissato continua
ad avere un costo $O(n)$. \\
Consideriamo, data una matrice $T_n$ in forma tridiagonale hermitiana, la seguente funzione
\[
\begin{array}{cccl}
 f: & \R  & \longto & \{ 0 , \ldots , n \}  \\
 & x & \longmapsto &  \sharp \{ \ \textrm{cambi di segno nel vettore} \ (p_0(x) , \ldots , p_n(x)) \ \}
\end{array}
\]
Perché $f$ sia ben definita conveniamo che se $p_j(x) = 0$ allora il suo segno è lo stesso di $p_{j-1}(x)$.
Questa è una buona definizione perché $p_0(x) = 1$ e quindi non è mai $0$.
\begin{os}
 Dalla definizione di $f$ si deduce che se il suo valore di cambia in un punto, allora in quel punto ci deve
 essere uno zero di (almeno) uno dei polinomi $p_j(x)$.
\end{os}
Osserviamo cosa succede quando un polinomio ``interno'', ovvero un $p_j$ con $j \in 1 \ldots n-1$ ha uno
zero in $\bar x$.
Innanzitutto sappiamo da quanto visto nel capitolo precedente che $p_{j-1}(x) \neq 0 \neq p_{j+1}(x)$; è
possibile però affinare questa affermazione osservando meglio la relazione ricorrente a tre termini:
\[
 p_{j+1}(x) = (\alpha_{j+1} - x) p_j(x) - \beta_j^2 p_{j-1}(x)
\]
Se $p_j(x) = 0$ si ottiene che $p_{j-1}(x) \cdot p_{j+1}(x) < 0$, e quindi anche in un intorno di $\bar x$.
Ricordando le definizione di $f$ si osserva che il segno di $p_j(x)$ è quindi indifferente ai fini del numero
di cambiamenti di segno. In particolare $f$ non cambia segno in $\bar x$.
Evidentemente $p_0(x)$ non può cambiare segno, quindi ci rimane da verificare cosa succede quando cambia segno
$p_n(x)$. \\
Ricordiamo alcune particolarità dei polinomi $p_n$ viste precedentemente
\begin{enumerate}[(i)]
 \item Ogni $p_j$ ha limite a $-\infty$ uguale a $+\infty$;
 \item Gli zeri di $p_{j-1}$ separano gli zeri di $p_j$, nel senso visto alla fine del precedente capitolo
 (Sezione~\ref{sec:eqsecolare});
 \item Tutti gli zeri di un polinomio di una matrice hermitiana tridiagonale sono semplici;
\end{enumerate}
Osserviamo quindi che per ogni radice di $p_n$ si ha che il segno della sua derivata e quello di $p_{n-1}$ sono
opposti. La condizione di semplicità sugli zeri ci assicura che la derivata sarà infatti diversa da $0$ in
$\bar x$ e quindi $p_n'(\bar x) p_{n-1}(\bar x) < 0$. In particolare la costanza locale del segno della derivata
intorno allo $0$ ci assicura anche la costanza di quello di $p_{j-1}$ e quindi il \emph{cambio di valore di $f$}.
In sintesi $f$ agisce come una sorta di ``contatore di zeri''. Dato infatti un qualsiasi intervallo $[a,b]$ si verifica,
usando le considerazioni precedenti che
\[
 f(b) - f(a) = \sharp \{ \ \text{zeri del polinomio caratteristico in } [a,b) \  \}
\]
 
Possiamo dunque applicare questa considerazione per applicare il metodo di bisezione per calcolare
uno specifico autovalore di $T_n$. Queste proprietà di $f$ ci permettono di calcolare infatti singolarmente
l'$i$-esimo autovalore di $T_n$ (ordinati con l'ordinamento di $\R$) senza bisogno di calcolare tutti gli altri.
Ricordando la diseguaglianza di Hirsch ($||\lambda|| \leq ||T_n||$) e supponendo di aver fissato $i$ possiamo
applicare il seguente procedimento.
 
Cominciamo col calcolare $f(0)$; il suo valore ci dirà quanti autovalori ci sono prima dello $0$, e quindi
nell'intervallo $[-||T_n||, 0)$. Se il nostro autovalore sta lì (ovvero $i \leq f(0)$) allora calcoliamo
$f(-\frac{||T_n||}{2})$ e ripetiamo il procedimento; in caso contrario procediamo analogamente con $f(\frac{||T_n||}{2})$. \\
Iterando un numero sufficiente di volte questo procedimento otterremo alla fine il valore desiderato.
 
Osserviamo che il procedimento può essere applicato a ad una qualsiasi matrice hermitiana, a patto
di portarla in forma tridiagonale (che è sempre possibile tramite matrici unitarie).
Ricordando le stime dei costi computazionali fatte in precedenza ci rendiamo conto che il costo principale è
quello della riduzione ($O(n^3)$) rispetto a quello del calcolo dell'autovalore ($O(n^2)$).
Se si volessero calcolare tutti gli autovalori il costo rimarrebbe ancora dell'ordine di $n^3$. Nonostante questo
esistono metodi (che vedremo in seguito) che permettono di calcolare tutti gli autovalori con una complessità di
$O(n^2)$ e che quindi saranno più interessanti per il calcolo dello spettro completo della matrice. \\
Il metodo di Sturm, comunque, trova delle applicazioni anche in questo ultimo caso nell'ambito del calcolo parallelo,
dove altri algoritmi non possono essere implementati. Si può osservare invece che, una volta ottenuta la fattorizzazione,
non sia difficile implementare il metodo di Sturm dividendo fra i vari processori gli autovalori da calcolare.
 
\section{Il metodo QR}
In questa sezione esporremo il metodo QR per il calcolo degli autovalori, che è il metodo più utilizzato
nelle moderne applicazioni del calcolo degli autovalori. Si basa sulla fattorizzazione QR già vista nel corso
di Analisi Numerica per la risoluzione di sistemi lineari. Questa può venire realizzata tramite matrici
di Householder. Cominceremo ricordando alcuni risultati riguardo la fattorizzazione, per poi esporre ed
analizzare il metodo e le sue implementazioni.
 
\subsection{La fattorizzazione QR}
\`E noto che ogni matrice $A \in \mat{\C}{n}$ si può fattorizzare nel seguente modo
\[
  A = QR
\]
dove $Q$ è una matrice unitaria e $R$ è una matrice triangolare superiore.
\begin{os}
 La fattorizzazione non è unica. Se supponiamo $A = QR$ e $S$ una matrice di fase, ovvero diagonale tale che $|s_{ii}|=1$
per ogni $i = 1 \ldots n$, allora $A = QS\herm{S}R$ è ancora un fattorizzazione. Se $S \neq I$ le fattorizzazioni sono
diverse. Si può però mostrare che non ci sono altre matrici (non di fase) per cui questo è vero. Si dice che
la fattorizzazione QR è \emph{essenzialmente unica}.
\end{os}
Il nostro scopo è costruire una successione di matrici tramite questo procedimento di fattorizzazione
che ci porti a determinare gli autovalori.
 
\subsection{Costruzione della successione}
Data una matrice $A$ di cui vogliamo conoscere gli autovalori, consideriamo
la successione definita nel seguente modo
\[
\left\{ \begin{array}{ll}
 A_0 & = A \\
 A_{k+1} & = R_k Q_k \quad \text{dove} \ Q_k R_k \ \text{è la fattorizzazione QR di} \ A_k
\end{array} \right.
\]
\begin{os} \label{os:qr_simili_ak}
 Osserviamo che per ogni $k$ $A_{k+1}$ è simile ad $A_k$. Infatti $A_{k+1} = \herm{Q_k}Q_k R_k Q_k$.
 Per ogni $k$ la matrice $A_k$ è simile a $A_{k+1}$ tramite trasformazione unitaria, che preserva
 il condizionamento del problema di calcolo degli autovalori. Questa quindi è una ``buona'' successione
 nel senso in cui ne avevamo parlato nella Sezione~\ref{sec:analisidelcondizionamento}.
\end{os}
 
Cominceremo ad analizzare il metodo QR facendo delle supposizioni piuttosto restrittive, che poi allenteremo
in seguito. Cominciamo col supporre che gli autovalori della matrice $A$ siano tutti distinti e ordinabili
per modulo in modo strettamente crescente, ovvero $|\lambda_1| < \ldots < |\lambda_n|$. Questo in
particolare implica che la matrice $A$ è diagonalizzabile e quindi esiste una $X$ invertibile tale che
$A = XDX^{-1}$. Supponiamo ora che $X^{-1}$ ammetta fattorizzazione $ X^{-1} = LU$\footnote{dove supponiamo
che $L$ sia triangolare inferiore con gli elementi delle diagonale uguali a $1$ e $U$ triangolare superiore}.
\begin{pr} \label{pr:metpot:ak}
 Se si ha la successione di matrici $A_k$ come quella definita sopra, e per ogni $k$ si considera $Q_k R_k$
 la\footnote{al solito, sarebbe corretto dire \textbf{una} fattorizzazione.} fattorizzazione QR di $A_k$, allora
 \[
  A^k = Q_0 Q_1 \ldots Q_{k-1} R_{k-1} \ldots R_1 R_0
 \]
\end{pr}
\begin{proof}
 Proviamo la tesi per induzione su $k$. Se $k = 1$ si ha $A^1 = A = Q_0 R_0$ che è banalmente vero. Se considero
 la tesi vera per $k$. Considerando le seguenti uguaglianze
 \[ \begin{array}{ll}
  \displaystyle
  A^{k+1} = \prod_{i=1}^{k+1} Q_0 R_0 = Q_0 (\prod_{i=1}^{k} R_0 Q_0) R_0 = Q_0 (\prod_{i=1}^{k} A_1) R_0 = \\
  \displaystyle
  = Q_0 (\prod_{i=1}^{k-1} Q_1 R_1) R_0 = Q_0 \prod_{i=1}^{k} Q_i \prod_{i=1}^{k} R_{k-i+1} R_0
 \end{array}
 \]
 si ottiene esattamente la tesi. Per quanto oscure possano sembrare provare a fare il calcolo con $k = 3$ o $4$
 potrebbe chiarire molto le idee.
\end{proof}
Ora vorremmo usare quanto scoperto sulle potenze di $A$ per studiare la convergenza del nostro metodo. Consideriamo
che $A^k = X D^{k} X^{-1}$ e quindi ricordando che esiste la fattorizzazione $LU$ di $X^{-1}$ si ha
\begin{equation}
A^k = XD^{k}LU= X D^k LD^{-k}D^k U
\end{equation}
Osserviamo ora che la matrice $D^k L D^{-k}$ è triangolare inferiore ed ha la diagonale con soli $1$. Se chiamiamo
$X = QR$ la fattorizzazione QR di $X$ possiamo scrivere
\begin{equation} \label{eq:metpot:1}
 A^k = X[I + \Gamma^{(k)}]D^k U = QR [ I + \Gamma^{(k)}] D^k U = Q[I + R\Gamma^{(k)}R^{-1}]RD^k U
\end{equation}
A questo punto consideriamo anche che per ogni $k$ esiste la fattorizzazione QR di $[I + R\Gamma^{(k)}R^{-1}] = P_k T_k$
ed in particolare possiamo scegliere $P_k$ a termini positivi. Osserviamo ora che i termini di $\Gamma^{(k)}$
sono dati dalla seguente relazione
\[
\gamma_{ij} =
 \left\{ \begin{array}{ll}
          0 & \text{se} \ i \leq j \\
          (\frac{\lambda_j}{\lambda_i})^{k} l_{ij} & \text{se} \ i > j
         \end{array} \right.
\]
ed in particolare ricordando che se $j < i$ si ha che $\lambda_j < \lambda_i$ si ottiene che la matrice $\Gamma^{(k)}$
tende ad una matrice diagonale per $k$ che tende all'infinito. Più precisamente, $\lim_{k \to \infty} \Gamma^{(k)} = I$.
Da questo si ottiene che $P_k \to I$ e anche $T_k \to I$\footnote{questo non sarebbe vero a priori, in quanto la scomposizione
QR è sempre definita a meno di una matrice di fase. Richiedere però che $P_k$ abbia elementi positivi ci permette
sia di definire univocamente la scomposizione desiderata sia di avere l'esistenza del limite.}
Confrontando l'equazione~\ref{eq:metpot:1} e la Proposizione~\ref{pr:metpot:ak} si ottengono due fattorizzazioni
QR della matrice $A^k$.
\[
 A^k = QP_k T_k R D^k U = Q_0 \ldots Q_k R_k \ldots R_0
\]
Due fattorizzazioni QR della stessa matrice devono forzatamente differire per una matrice di fase e quindi
si ottengono le due relazioni
\[
 \left\{ \begin{array}{l}
          Q_0 \ldots Q_{k-1} = Q P_k S_k \\
	  R_{k-1} \ldots R_0 = \herm{S_k}T_k R D^k U
         \end{array} \right.
\]
Possiamo riscrivere ora $Q_k$ ed $R_k$ in modo da riuscire ad utilizzare queste relazioni\footnote{Si può
osservare quanto i passaggi di questa dimostrazione siano la cosa meno intuitiva pensabile (o quasi). Presumibilmente
questa è risultato di anni di affinamento e di ``pulizia'' }
\[
 \left\{ \begin{array}{l}
         Q_k = \herm{(Q_0 \ldots Q_{k-1})}(Q_0 \ldots Q_k) = \herm{S_{k}}\herm{P_{k}}\herm{Q}Q P_{k+1} S_{k+1} =
         \herm{S_{k}}\herm{P_{k}}P_{k+1} S_{k+1}\\
         R_k = (R_k \ldots R_0)(R_{k-1} \ldots R_0)^{-1} = \herm{S_{k+1}}T_{k+1} R D^{k+1} U U^{-1} D^{-k} R^{-1} T_{k}^{-1} S_{k}
        \end{array} \right.
\]
Dunque $Q_k R_k = A_k = \herm{S_{k}}\herm{P_{k}}P_{k+1} T_{k+1} R D R^{-1} T_{k}^{-1} S_{k}$ e ricordando che
$T_k$ e $P_k$ al limite vanno all'identità se scriviamo la nostra uguaglianza per $k~\to~\infty$ otteniamo
$S_{k} A_k \herm{S_{k}} = R D R^{-1}$. Possiamo quindi osservare che gli elementi sulla diagonale di $RDR^{-1}$
sono gli stessi di $D$ (grazie al fatto che $R$ è triangolare superiore). In particolare gli elementi diagonali
sono gli autovalori di $A$ e quindi abbiamo provato che il metodo converge.
\begin{os}
 In realtà non è rilevante conoscere esplicitamente $S_k$ perché siamo interessati solo a conoscere gli elementi
 sulla diagonale di $RDR^{-1}$. Avendo che gli elementi di $S_k$ sono di modulo $1$ gli elementi diagonali vengono
 moltiplicati per un numero complesso e per il suo coniugato, lasciandoli invariati.
\end{os}
 
\subsection{Il costo computazionale} \label{subsec:qr_costo}
Vorremmo ora valutare il costo computazionale di questo metodo. Consideriamo dapprima il caso di una matrice
generale. Sappiamo che il costo di una fattorizzazione $QR$ è $O(n^3)$ e sperimentalmente si ottiene che il
numero di passi per avere convergenza cresce linearmente con la dimensione. Questo ci permette di concludere
che in un caso totalmente generale il costo computazionale del metodo QR è $O(n^4)$.
 
Facciamo ora qualche supposizione sulla matrice. Prendiamo ad esempio una matrice tridiagonale hermitiana.
In questo caso ci ricordiamo che il costo del calcolo della fattorizzazione QR (almeno della prima) è $O(n)$.
Possiamo osservare che la matrice ottenuta dopo il primo passo è ancora tridiagonale hermitiana, e quindi
le considerazioni fatte sul primo passo valgono anche per i seguenti. In particolare il costo del metodo
è $O(n^2)$.
 
Analogamente per le matrici di Hessenberg (anche se non lo mostriamo ora) il costo è di $O(n^3)$ perché
il calcolo della scomposizione QR è $O(n^2)$.
 
\subsection{Migliorare la convergenza}
Il metodo QR come esposto fino ad ora funziona ma ha il difetto di avere una convergenza di tipo
lineare. In generale questo può essere un problema perché più gli autovalori della matrice sono
vicini più, potenzialmente, potrebbe essere lenta (o addirittura assente) la convergenza.
 
Possiamo osservare meglio questo con un esempio; consideriamo una matrice $2\times 2$ con un
elemento ``piccolo'' in posizione $(2,1)$, ovvero sulla buona strada per essere triangolarizzata\footnote{%
Ovviamente tutta questa frase è terribilmente imprecisa, ma il nostro scopo è, per ora, solo dare
un'idea dei vantaggi che si potrebbero ottenere da un nuovo approccia e non formalizzarlo in
alcun modo}.
Proviamo a vedere cosa succede effettuando un passo generico dell'iterazione. Consideriamo la matrice $A_k$
con un elemento piccolo $\eps$ sotto la diagonale e una sua scomposizione QR
\[
 A_k = \left[ \begin{array}{cc}
                a & b \\
                \eps & c \\
               \end{array}\right]
 = \left[ \begin{array}{cc}
           \frac{a}{\sqrt{a^2 + \eps^2}} & \frac{-\eps}{\sqrt{a^2+\eps^2}} \\
           \frac{\eps}{\sqrt{a^2+\eps^2}} & \frac{a}{\sqrt{a^2 + \eps^2}}
          \end{array} \right]
  \left[ \begin{array}{cc}
          \sqrt{a^2 + \eps^2} & \frac{ab + \eps c}{\sqrt{a^2 + \eps^2}} \\
	  0 & \frac{ac - \eps b}{\sqrt{a^2+\eps^2}}
         \end{array} \right]  = Q_kR_k
\]
Si avrà che
\[
 A_{k+1} = R_kQ_k =   \left[ \begin{array}{cc}
          \sqrt{a^2 + \eps^2} & \frac{ab + \eps c}{{a^2 + \eps^2}} \\
	  0 & \frac{ac - \eps b}{{a^2+\eps^2}}
         \end{array} \right]
 \left[ \begin{array}{cc}
           \frac{a}{\sqrt{a^2 + \eps^2}} & \frac{-\eps}{\sqrt{a^2+\eps^2}} \\
           \frac{\eps}{\sqrt{a^2+\eps^2}} & \frac{a}{\sqrt{a^2 + \eps^2}}
          \end{array} \right]
= \left[ \begin{array}{cc}
          \times & \times \\
	  \frac{\eps(ac - \eps b)}{a^2 + \eps^2} & \times
         \end{array} \right]
\]
Considerando che $ac -\eps b$ è il determinante delle matrice (e quindi
il prodotto degli autovalori) e $a^2 + \eps^2$  una buona
approssimazione del quadrato di un autovalore\footnote{questo perché il vettore $e_1$ va a finire in se stesso moltiplicato
per $a$, a meno di epsilon che è piccolo}, $\eps$ è stato diminuito di un fattore dell'ordine
di $\frac{\lambda_1}{\lambda_2}$, come ci si aspettava.
 
Consideriamo ora cosa succede applicando una tecnica di \emph{shifting}. Calcoliamo la fattorizzazione QR della
matrice $A-cI$, invertiamo l'ordine dei fattori e risommiamo $cI$.
Abbiamo dunque
\[
  A_k - cI = \hat Q \hat R \quad \longrightarrow \quad A_{k+1} = \hat R \hat Q + cI
\]
\begin{os}
 $A_k$ ed $A_{k+1}$ sono simili. Consideriamo infatti le seguenti uguaglianze
 \[
  A_{k+1} = \hat R \hat Q + cI = \herm{\hat Q}\hat Q \hat R \hat Q + cI = \herm{\hat Q} (A_k - cI) \hat Q + cI
  = \herm{\hat Q} A_k \hat Q
 \]
 E quindi le due matrici sono simili per trasformazione hermitiana, il che preserva anche il condizionamento del
 problema come visto nella Sezione~\ref{sec:analisidelcondizionamento}.
\end{os}
 
Possiamo osservare che risvolgendo i conti con il ``nuovo'' algoritmo si ottiene in posizione sottodiagonale
un fattore dell'ordine di $\eps^2$. Sembra quindi che questo nuovo metodo abbia una convergenza molto più veloce
del precedente!
 
Anche se questo esempio ci fornisce solo delle osservazioni in un caso particolare e non una
dimostrazioni formale delle proprietà di convergenza che si potrebbero ottenere in generale
con uno shifting, possiamo in ogni caso discutere come applicarlo. Esistono teoremi che provano
l'effettiva efficienza di questo approccio ma sono piuttosto tecnici e ``conterecci'' (cit.) e
non li affronteremo in questo corso.
\begin{figure}[hb]
\[
 A_k = \left[ \begin{array}{cccc}
               \alpha_{11}^{(k)} & \times & \hdots & \times \\
	       \beta_{1}^k & \ddots & & \vdots \\
                           & \ddots & \ddots& \times \\
                & & \beta_{n-1}^k & \alpha_{nn}^{(k)}
              \end{array} \right]
\]
\caption{Struttura di una matrice in forma di Hessenberg superiore}
\label{fig:hessenbergsup}
\end{figure}
Appurato che questo trucco funziona, ci chiediamo come applicarlo ad una matrice generale. L'idea
è di togliere alla matrice $A_k$ un'altra matrice $\alpha_k I$, dove $\alpha_k = \alpha_{nn}^{(k)}$
(si veda la Figura~\ref{fig:hessenbergsup}).
Per semplicità osserviamo il caso delle matrici di Hessenberg, che poi sarà il caso implementato in pratica
date le considerazioni sulla complessità fatte in precedenza.
 
Il caso di una matrice complessa viene gestito esattamente con lo shifting di $-\alpha_{nn}^{(k)}I$ come descritto
sopra, e si verifica che le cose funzionano bene. Più complesso, invece, è il caso di una matrice reale.
In quella situazione infatti gli autovalori possono essere complessi coniugati e non si desidera passare da aritmetica
reale ad aritmetica complessa. Questo porterebbe a problemi di complessità, ma il danno più grave sarebbe
perdere la simmetria\footnote{Per simmetria intendiamo la certezza di ottenere autovalori complessi coniugati,
cosa di cui non si potrebbe essere certi con una matrice complessa qualsiasi. } della matrice,
e quindi in definitiva parecchia informazione sul problema.
\begin{os}
 Nel caso complesso Hessenberg risulta semplice determinare la qualità dell'approssimazione ottenuta
 dell'autovalore. Possiamo assumere di fermare l'iterazione quando
 \[
  |\beta_n| \leq u( |\alpha_{nn}^{(k)}| + |\alpha_{n-1,n-1}^{(k)}|)
 \]
 dove $u$ è la precisione di macchina,
 ed a quel punto assumere che $\alpha_{nn}^{(k)}$ è una buona approssimazione dell'autovalore
 cercato e ridurci a ripetere il procedimento sul minore di nord ovest di ordine $n-1$, che è ancora in forma
 di Hessenberg.
\end{os}
Per determinare come eseguire lo shifting nel caso reale conviene cambiare punto di vista su quanto effettuato fino
ad ora. Possiamo osservare che se poniamo $p_k(z) = z - \alpha_k$ di fatto lo shifting è equivalente a considerare
la matrice $\hat A_k = p(A_k)$. Dato che riteniamo che $\alpha_k$ sia una buona approssimazione dell'autovalore
cercato della matrice per analogia nel caso reale possiamo utilizzare $p_k(x) = (x - \lambda)(x - \con{\lambda})$
che è un polinomio reale. Si verifica (e, ancora una volta, noi non lo faremo) che tutto continua a funzionare
bene ma sorgono dei problemi di complessità. \\
Il polinomio $p_k$ è infatti adesso di II grado e non è per nulla banale, in generale, calcolare $A_k^2$.
Inoltre $p_k(A_k)$ non è neppure una matrice in forma di Hessenberg, il che sembrerebbe farci perdere tutti
i vantaggi computazionali che avevamo fino a questo momento. Per fortuna esiste una tecnica che, modificando
l'iterazione, riesce a riolvere questo problema.
 
\subsection{La tecnica del Bulge Chasing}
Nell'applicazione della tecnica dello shifting vista nel paragrafo precedente, abbiamo supposto di
calcolare $p_k(A_k)$ e dopo calcolare la fattorizzazione QR. Ci siamo però accorto che se $p_k$ è
un polinomio di secondo grado vengono persi tutti i vantaggi computazionali dati dalla forma di Hessenberg.
Vogliamo quindi ora modificare la fattorizzazione QR \textbf{mentre} calcoliamo $p_k(A_k)$ in modo da
rendere la risultante matrice in forma di Hessenberg e diminuire la complessità del calcolo del polinomio.
La fattorizzazione QR ottenuta sarà ``essenzialmente uguale''\footnote{ovvero sarà uguale a meno di una
matrice di fase reale} a quella che avremmo ottenuto seguendo
il procedimento precedente; questo modo di operare è detto \emph{Bulge Chasing}\footnote{bulge chasing in inglese significa inseguimento
del bernoccolo, e s riesce ad apprezzare il senso del nome osservando come questa tecnica si traduce graficamente sugli elementi
non zero della matrice. Questi si comportano infatti come un bernoccolo che scende sempre di più sulla
sottodiagonale per poi sparire lasciando una matrice in forma di Hessenberg superiore.}.
 
Possiamo cominciare a calcolare solo la prima colonna della matrice $p_k(A_k)$, ovvero $p_k(A_k)e_1$.
Una volta calcolata questa possiamo calcolare anche una matrice di Householder $P_0$ tale che
$P_0 p_k(A_K)e_1 = \alpha e_1$. \\
Si costruiscono poi altre $n-1$ matrici di Householder tali che fissato
\[
 Z = P_0 \ldots P_{n-1}
\]
si abbia
\[
 A_{k+1} = \herm{Z} A_k Z
\]
ed è possibile dimostrare che questa iterazione produce una matrice simile a meno di una matrice di fase reale
rispetto a quella presentata in precedenza.
 
\section{Divide et Impera}
\subsection{L'idea}
Introdurremo ora un altro metodo, piuttosto recente, per il calcolo di tutti gli autovalori della matrice.
Questo è stato introdotto nel 1980 da Cuppen. Il metodo permette il calcolo di autovalori solo per matrici
tridiagonali hermitiane e l'osservazione che ne sta alla base è che una matrice tridiagonale è ``quasi''
diagonale a blocchi, più precisamente se $H$ è la nostra matrice
\[
 H = \left[ \begin{array}{cc|cc}
             \sblocko{T_1}{2} & &  \\
            & & & \\ \hline
             & & \sblocke{T_2}{2} \\
            & & & \\
            \end{array} \right]
+ \left[ \begin{array}{cc|cc}
          & & & \\
          & \beta_{n-1} & & \\ \hline
          & & \beta_{n-1} & \\
          & & &
         \end{array} \right]
\]
Ovvero è una diagonale a blocchi più una correzione di rango 2. A ben vedere, però, si può fare anche di meglio.
Se sostituiamo in $T_1$ il valore $\alpha_{\frac{n}{2}}$ con $\hat \alpha_{\frac{n}{2}} = \alpha_{\frac{n}{2}} - \beta_{n-1}$
e in $T_2$ il valore $\alpha_{\frac{n}{2} + 1}$ con $\hat \alpha_{\frac{n}{2} + 1} = \alpha_{\frac n 2 + 1} - \beta_{n-1}$, otteniamo
una matrice a blocchi tridiagonali più un correzione di rango $1$:
\[
 H = \left[ \begin{array}{cc|cc}
             \sblocko{\hat T_1}{2} & &  \\
            & & & \\ \hline
             & & \sblocke{\hat T_2}{2} \\
            & & & \\
            \end{array} \right]
+ \left[ \begin{array}{cc|cc}
          & & & \\
          & \beta_{n-1} & \beta_{n-1}& \\ \hline
          & \beta_{n-1} & \beta_{n-1} & \\
          & & &
         \end{array} \right]
\]
Supponiamo ora di conoscere gli autovalori di $\hat T_1$ e di $\hat T_2$, ed in particolare le matrici
unitarie $Q_1$ e $Q_2$ che ci permettono di diagonalizzarle, ovvero
\[
 T_1 = Q_1 D_1 \herm Q_1 \qquad
 T_2 = Q_2 D_2 \herm Q_2
\]
con $D_1$ e $D_2$ le matrici con gli autovalori sulla diagonale e $0$ altrove.
 
Tutto questo può funzionare se
riusciamo a mostrare che siamo veramente in grado di calcolare facilmente gli autovalori della matrice a blocchi
con la correzione.
 
Se consideriamo ora la matrice $D$ così formata
\[
 \hat D = \left[ \begin{array}{cc|cc}
                  \sblocko{D_1}{2} & & \\
                  & & & \\ \hline
                  & & \sblocke{D_2}{2} \\
                  & & &
                 \end{array} \right]
\]
e calcoliamo $B = \herm Q A Q$ otteniamo
\[
 B = \herm Q A Q = \hat D + \beta_{\frac{n}{2}}z\herm z
\]
con $z = \left[ \begin{array}{c} q_1 \\ q_2 \end{array}\right]$ dove $q_i$ è l'ultima colonna di $\herm Q_i$
per $i = 1,2$.
 
\subsection{La ricerca degli autovalori}
Abbiamo quindi ricondotto il problema a calcolare gli autovalori di una matrice $D + w\herm w$ dove $D$
è diagonale e $w$ è un vettore qualsiasi. Calcolando il polinomio caratteristico di $B$ ed assumendo
che $\lambda \neq \hat d_i$  otteniamo
\[
  p(\lambda) = \deter{\lambda I - B} = \deter{ \lambda I - D }\cdot \deter{I + \theta(\lambda I - D)^{-1} z\herm z}
\]
Osserviamo che stiamo assumendo che $\lambda I - D$ sia invertibile, e quindi questa formula per il calcolo
del polinomio caratteristico vale solo per $\lambda \in \C \setminus (\spe{D_1} \cup \spe{D_2})$. Sapendo che $p$
è continuo possiamo però concludere che vale per tutto $\C$. \\
Notiamo inoltre che $( I + \theta(\lambda I - D)^{-1}z\herm z)$ è una matrice elementare e quindi ha determinante\footnote{%
Ricordiamo che in genearale una matrice elementare $I - \sigma u \herm v$ ha determinante $1 - \sigma \herm v u$}
$1 - \theta (\lambda I - D)^{-1} \herm zz$; sviluppando ulteriormente si ottiene
\begin{equation} \label{eqn:polmindivetimp}
 p(\lambda) = \prod_{j=1}^{n} (\lambda - \hat d_j) ( 1 + \theta (\lambda I - D)^{-1} \herm zz )
 = \prod_{j=1}^{n} (\lambda - \hat d_j) ( 1 + \theta \sum_{k=1}^{n} \frac{z_k^2}{\lambda - \hat d_k} )
\end{equation}
ed estendendo per continuità a tutto $\C$ si ottiene infine
\[
 p(\lambda) = \prod_{j=1}^n (\lambda - \hat d_j) + \theta\sum_{k=1}^n z_k \prod_{\substack{s=1 \\ s\neq k}}^n (\lambda - \hat d_s)
\]
che è una sorta di equazione secolare come quella che si era ottenuta partendo dalle relazioni ricorrenti a tre
termini nella Sezione~\ref{sec:eqsecolare}.
 
Siamo ora costretti a fare delle assunzioni sulla struttura di $z$ e di $\hat D$. In realtà ci renderemo conto
in seguito che queste non sono restrittive perché nei casi in cui non saranno verificate ci troveremo con un
problema semplificato, che potremo riportare a questo caso.\footnote{Penso, perché in realtà non mi è ancora
molto chiaro come si dovrebbe fare}. \\
Assumiamo dunque che
\begin{itemize}
 \item $z_j \neq 0$ per ogni $j = 1 \ldots n$;
 \item $\hat d_k \neq \hat d_j$ per ogni $k \neq j$, ovvero che tutti gli autovalori siano distinti;
\end{itemize}
A questo punto dalla relazione~\ref{eqn:polmindivetimp} ricaviamo che
\begin{equation} \label{eqn:eqsecdivetimp}
 p(\lambda) = 0 \iff 1 + \theta \sum_{j=1}^n \frac{z_j^2}{\lambda - \hat d_j} = 0
\end{equation}
 
e quindi gli autovalori delle matrici di ordine $\frac{n}{2}$ separano quelli della matrice grande, e si
può applicare (ad esempio) il metodo di bisezione, o in generale altri procedimenti di iterazione
funzionale.
 
\begin{os}
 Notiamo che per applicare la bisezione nella (\ref{eqn:eqsecdivetimp}) abbiamo bisogno di conoscere, oltre
 ai $\hat d_j$ anche gli $z_j$ e ricordando che $z$ è costruito a partire dalle matrici $Q_1$ e $Q_2$, è chiaro
 che dobbiamo conoscere anche queste due. Queste erano le matrici che diagonalizzavano $T_1$ e $T_2$ e si
 ottenevano quindi \textbf{conoscendo gli autovettori} di $T_1$ e $T_2$. Dobbiamo allora, per poter procedere
 con il metodo, calcolarci tutti gli autovettori.
\end{os}
 
\subsection{Autovettori e problemi di stabilità}
Il fatto di dover calcolare gli autovettori ha un lato positivo ed uno negativo, precisamente
\begin{enumerate}[(a)]
 \item Il calcolo effettivo risulta inaspettatamente semplice, tanto che, trovati gli autovalori saremo
 capaci di trovare ogni autovettore con costo $O(n)$ e quindi di ottenerli tutti con costo $O(n^2)$;
 \item Il problema del calcolo degli autovettori non è in generale ben condizionato e non si riescono a
 dare condizioni perchè lo sia. Dato che influenzerà anche il calcolo degli autovalori nei passaggi successivi,
 non possiamo garantire l'accuratezza del risultato finale;
\end{enumerate}
 
L'idea per trovare gli autovettori è cercare quelli di $\hat D + \beta_{\frac n 2}z\herm z$ e poi una
volta ottenuta la matrice $U$ con questi, moltiplicarla per la matrice a blocchi con $Q_1$ e $Q_2$ sulla
diagonale (che d'ora in poi chiameremo $\hat Q$) per trovare gli autovettori della matrice $T$. \\
Osserviamo che se $v$ è autovettore per la matrice $D + \beta_{\frac n 2}z\herm z$ allora
\[
 (\hat D + \theta z \herm z - \lambda I)v = 0
\]
e quindi
\[
 (\hat D - \lambda I)v = -\theta z \herm z v
\]
Se fosse $\herm z v = 0$ avrei che $(\hat D - \lambda I)v = 0$ ma questo non è possibile perché $v$
è non nullo e $\hat D - \lambda I$ è invertibile (come assunto precedentemente). Quindi deve essere
$\herm z v \neq 0$ e quindi dato che $v$ è determinato a meno di uno scalare possiamo assumere
$\herm z v = 1$. In definitiva si ottiene
\[
 v = - \theta z (\hat D - \lambda I)^{-1}
\]
che si può scrivere anche più esplicitamente come
\begin{equation}  \label{eqn:calcoloautovettdivetimp}
 v_j = \frac{- \theta z}{d_j - \lambda}
\end{equation}
Riassumendo, vogliamo trovare la matrice $U$ con i $v_j$ sulla diagonale e questo si può fare con
la Formula~\ref{eqn:calcoloautovettdivetimp} con costo $O(n^2)$. Ottenuta la matrice $U$ avremo
la matrice per diagonalizzare $T$ calcolando $\hat Q U$, e potremmo quindi procedere con il metodo.
 
Sorgono però dei problemi computazionali non banali, più precisamente
\begin{enumerate}[(i)]
 \item La matrice $U$ deve essere ortogonale ma questo non è garantito dall'algoritmo. Se ad un passo
 abbiamo perdita di ortogonalità questo poi causerà un accumulo di errore in tutti i passi successivi
 portando ad una imprecisione del calcolo. Sembra che in una implementazione del 2006 si sia riusciti
 a risolvere questo problema. %% TODO: inserire una citazione
 \item Come già sottolineato prima il calcolo degli autovalori potrebbe essere mal condizionato dal principio
 portando ad un accumulo di errore notevole ad ogni passo.
 \item Apparentemente il costo computazionale del calcolo di $\hat Q U$ non è trascurabile in quanto
 un prodotto matrice per matrice in generale costa $O(n^3)$. Fortunatamente $U$ non è una matrice qualsiasi:
 matrici del tipo $u_{ij} = \frac{r_i s_j}{x_i - y_j}$ sono dette matrici \emph{Cauchy-like} ed esistono
 algoritmi per calcolare il prodotto fra una Cauchy-like ed una matrice generica con costo $O(n^2\log{n})$
\end{enumerate}
 
Come ultima nota, ricordiamo che questo algoritmo è applicabile solo alle matrici tridiagonali hermitiante
(e, previa tridagonalizzazione, a qualsiasi matrice hermitiana) ma non alle matrici in forma di Hessenberg
che non presentano purtroppo alcuna proprietà di separazione\footnote{ovvero quella che ci permette di trovare
i limite su cui applicare la bisezione o qualsiasi altra iterazione funzionale}. \\
Sono stati fatti molti tentativi a questo proposito, ma non esiste attualmente nessuna implementazione
robusta del metodo Divide et Impera sulle matrici di Hessenberg.
 
\section{Il metodo delle potenze} \label{sec:metodopotenze}
Per ultimo analizzeremo il metodo delle potenze che ci permette di trovare gli autovettori di una matrice
conoscendo gli autovalori. In generale è particolarmente comodo per calcolare l'autovalore dominante.
Consideriamo dunque questo caso per primo
 
\subsection{Il metodo diretto}
Supponiamo $A \in \mat{\C}{n}$ diagonalizzabile, per semplicità. Elimineremo questa ipotesi in seguito.
Consideriamo gli autovalori ordinati in modo che
\[
 |\lambda_1| \geq |\lambda_2| \geq \ldots \geq |\lambda_n|
\]
e aggiungiamo l'ipotesi che $\lambda_1$ sia semplice e che sia in modulo strettamente maggiore degli altri.
Il metodo si basa su questa semplice osservazione (la cui paternità viene, ancora una volta, attribuita a Gauss):
Se prendiamo un vettore $y_0 \in \C^n$ e la nostra scelta non è incredibilmente sfortunata\footnote{%
dove con incredibilmente sfortunata intendiamo che il prodotto scalare tra $y_0$ e l'autovettore relativo
a $\lambda_1$ è esattamente uguale a $0$} possiamo considerare la successione
\begin{equation}
 \left\{ \begin{array}{ll}
          y_0 & = y_0 \\
          y_k & = Ay_{k-1}
         \end{array} \right.
\end{equation}
Se scriviamo $y_0$ nella base di autovettori per $A$ abbiamo
\begin{equation}
 y_0 = \sum_{i=1}^{n} \alpha_i v_i
\end{equation}
e quindi osservando che $y_k = Ay_{k-1} = A^{k}y_0$ dopo $k$ iterazioni otteniamo
\begin{equation} \label{eqn:metpotite}
 y_k = A^ky_0 = \sum_{i=1}^n \alpha_i \lambda_i^k v_i =
 \lambda_i^k \cdot \Big( \alpha_1 v_1 + \underbrace{\sum_{i=2}^n \alpha_i (\frac{\lambda_i}{\lambda_1})^k v_i}_{%
\text{tende a}\ 0\ \text{se}\ k \to \infty} \Big)
\end{equation}
e ricordando l'ipotesi fatta precedentemente, ovvero che $|\lambda_1| > |\lambda_i|$ per ogni $i$ abbiamo
che la seconda parte dell'uguaglianza va a $0$ per $k \to \infty$ e quindi la successione di vettori
tende ad un multiplo di $v_1$, ovvero l'autovettore cercato.
 
Purtroppo questa implementazione del metodo non è realistica, perché si incorrerebbe sicuramente nel caso di
overflow o underflow in poche iterazioni, e probabilmente prima di avere una buona approssimazione dell'autovettore.
Si ricorre quindi al semplice stratagemma di riscalare ad ogni iterazione il vettore ottenuto, ad esempio
dividendolo per la sua norma infinito. In questo modo otteniamo una successione di vettori unitari che convergeranno
all' autovettore unitario. \\
Possiamo ora emprimere la successione in questo modo
\[
 \left\{ \begin{array}{ll}
          y_0 & = y_0 \\
	  y_k & = \displaystyle \frac{Ay_{k-1}}{||Ay_{k-1}||_\infty}
         \end{array} \right.
\]
\begin{os}
 Questo metodo mi permette anche di valutare e raffinare la stima dell'autovalore dominante. Se considero
 infatti il rapporto fra le $j$-esime componenti di due vettori consecutivi, ovvero $\frac{y_k,j}{y_{k-1,j}}$
 non è difficile verificare tramite la (\ref{eqn:metpotite}) che
 converge all'autovalore cercato.
\end{os}
 
Esistono delle varianti del metodo delle potenze che sfruttano le idee qui esposte per risolvere problemi simili
e migliorare le condizioni di convergenza.
 
\subsection{Il metodo inverso con shift}
Osserviamo come prima cosa che con pochi cambiamenti si potrebbe calcolare l'autovalore di modulo più piccolo
considerando, invece che la matrice $A$, la matrice $A^{-1}$. \\
Potremmo riscrivere l'iterazione in questo modo
\begin{equation}
 \left\{ \begin{array}{ll}
          y_0 & = y_0 \\
          y_k & = \beta_{k} \cdot A^{-1}y_{k-1}
         \end{array} \right.
\end{equation}
dove $\beta_{k} = ||y_{k}||_\infty$\footnote{adotteremo questa notazione per tutta la sezione, d'ora in poi}.
Questo procedimento ci imporrebbe però di calcolare l'inversa di $A$, un procedimento che in generale si tende ad
evitare per problemi di instabilità. Possiamo quindi rileggere l'espressione precedente come risoluzione di
un sistema lineare, ovvero
\begin{equation}
 \left\{ \begin{array}{lll}
          y_0 &=& y_0 \\
          Az_k &=& y_{k-1} \\
          y_k &=& \beta_k \cdot z_k
         \end{array} \right.
\end{equation}
A questo punto possiamo scegliere varie tecniche per risolvere il sistema lineare; sottolineamo che un metodo diretto
che utilizza una fattorizzazione in questo caso è molto conveniente rispetto ad uno iterativo. Se calcoliamo
la fattorizzazione di $A$, infatti, possiamo poi usarla per risolvere tutti i sistemi lineari della successione con
un costo basso.
 
In particolare, ricordiamo che per una tridiagonale hermitiana e per una matrice in forma di Hessenberg superiore
il costo di una fattorizzazione QR sono rispettivamente $O(n)$ e $O(n^2)$, e sono anche il costo della risoluzione
del sistema lineare. Si verifica sperimentalmente che, partendo da una buona approssimazione degli autovalori,
il metodo converge in genere in 3-4 passi, indipendemente dalla dimensione di $A$. Si conclude quindi che il calcolo
dell'autovettore dominante risulta essere $O(n)$ per le tridiagonali e $O(n^2)$ per le matrici in forma di
Hessenberg.
 
Ci poniamo ora il problema di come calcolare gli autovettori che stanno nella parte centrale dello spettro, ed è qui
che interviene lo \emph{shifting}. Non è difficile osservare che gli autovalori della matrice $A - \gamma I$ sono gli
stessi di $A$ ``shiftati'' di $\gamma$, ovvero sono
\[
 \spe{A - \gamma I} = \{ \lambda - \gamma \ | \ \lambda \in \spe{A} \ \}
\]
Supponiamo ora di avere una buona approssimazione $\hat \lambda_j$ dell'autovalore $\lambda_j$, dove con ``buona
approssimazione'' intendiamo tale che
\[
 |\hat \lambda_j - \lambda_j| < |\hat \lambda_j - \lambda_k| \quad \forall k \in \spe{A}\setminus\{\lambda_j\}
\]
In questo caso $\hat \lambda_j - \lambda$ è l'autovalore di modulo più piccolo della matrice $A - \hat \lambda_j I$
e il suo autovettore è quello relativo all'autovalore $\lambda_j$ della matrice $A$. Possiamo quindi applicare
il metodo delle potenze inverso a $A - \lambda_j I$ e ottenere l'autovettore cercato. Questa operazione prende
il nome di \emph{metodo delle potenze inverso con shift}.
 
\begin{os}
 In questo modo si può applicare un metodo qualsiasi per ottenere un'approssimazione degli autovalori (ad esempio
 il metodo QR) e poi calcolare tutti gli autovettori della matrice con il metodo delle potenze inverso con shift.
 Se la matrice è tridiagonale hermitiana si possono ottenere tutti con costo $O(n^2)$, mentre se è in forma di
 Hessenberg superiore il costo diventa $O(n^3)$ per le considerazioni fatte in precedenza.
\end{os}
 
\begin{os}
 Un'altra cosa che si può notare è che il metodo fornisce anche un possibile raffinamento per l'autovalore
 $\lambda_j$ tramite il rapporto fra le componenti di $y_k$ e di $y_{k-1}$. Si potrebbe essere tentati
 di sostituire questo eventuale raffinamento al posto dello shifting ma si verifica che questa non è necessariamente
 una buona pratica. Esistono alcuni casi in cui può inficiare la convergenza del metodo.
\end{os}