\chapter{Risoluzione di sistemi lineari}
 
In questo capitolo affronteremo il tema della risoluzione dei sistemi lineari introducendo dei metodi
iterativi specifici per alcune classi di matrici particolarmente interessanti; queste individueranno
solitamente sistemi troppo grandi per essere risolti tramite metodi diretti e in cui la convergenza
dei metodo classici (Jacobi e Gauss Seidel) è quasi assente.
 
\section{Sitemi lineari con matrici definite positive}
 
\subsection{Metodo del gradiente}
Supponiamo di avere una matrice $A \in \mat{\R}{n}$ definita positiva ed un sistema lineare
\[
 Ax = b
\]
In molti casi pratici (come ad esempio nello studio delle Vibrazioni, vedi Sezione~\ref{sec:vibsistcont})
ci si trova a risolvere sistemi lineari con matrici definite positive sparse o strutturate\footnote{in
cui in generale il costo del prodotto matrice vettore è basso} e in cui è conveniente usare un metodo
iterativo.
 
Consideriamo la seguente funzione
\begin{equation}
 \Phi(x) = \frac{1}{2} \trasp x A x - \trasp x b
\end{equation}
Osserviamo che il suo gradiente è dato dalla seguente espressione
\begin{equation}
 \nabla \Phi(x) = Ax - b
\end{equation}
e quindi vale $0$ solo se $x$ è una soluzione del sistema lineare. Possiamo quindi riformulare la ricerca
della soluzione del sistema lineare come ricerca dei punti stazionari della funzione $\Phi$.
Possiamo osservare che se $\bar x$ è un punto stazionario per $\Phi$ allora calcolando la matrice
delle derivate seconde si ottiene esattamente $A$; ricordando che $A$ è definita positiva si può concludere
che $\bar x$ è un punto di minimo per $\Phi$ e quindi
\[
 A\bar x = b \iff \bar x = \min_{x \in \R^n} ( \Phi(x) )
\]
In sintesi abbiamo ricondotto il problema della risoluzione del sistema lineare ad un problema di minimizzazione.
Ricordando che siamo alla ricerca di un metodo iterativo vorremmo affrontare il problema nel modo seguente
\begin{enumerate}[(a)]
 \item Scegliamo un punto $x_0$ a caso;
 \item Cerchiamo di determinare in che direzione si trova il minimo;
 \item Ci spostiamo ponendo $x_1 = x_0 + \alpha_0v_0$ dove $v_0$ è la \emph{direzione di decrescita}
 appena determinata e $\alpha_0$ un opportuno scalare;
\end{enumerate}
 
I \emph{metodi del gradiente}, ovvero quei metodi basati sulle considerazioni appena fatte, assumono
nomi diversi a seconda della tecnica scelta per determinare la direzione di decrescita.
 
\subsection{Il metodo del gradiente ottimo}
Una prima scelta piuttosto naturale per la direzione di decrescita nel punto $x_k$ potrebbe essere $-\nabla \Phi(x_k)$.
Ricordiamo infatti dall'analisi che il gradiente indica la direzione in cui la funzione ``cresce'' di più.
Questa scelta si dice del \emph{gradiente ottimo} ed è stata, storicamente, la prima ad essere implementata
e studiata. \\
Una volta scelta la direzione dobbiamo determinare $\alpha_k$. Per fare questo studiamo la seguente funzione di
$\alpha$ che valuta $\Phi$ sulla retta $x_k + \alpha v_k$, con $\alpha \in \R$:
\[
 g(\alpha) = \Phi(x_k + \alpha v_k)
\]
Ricordando che vogliamo trovare il minimo di $\Phi$ e osservando che $g$ è convessa cerchiamo anche qui un punto
stazionario di $g$; questo sarà l'$\alpha_k$ che ci permette di ottenere il valore minimo di $\Phi$ sulla direzione
determinata da $v_k$. \\
Otteniamo
\[
 g'(\alpha) = \alpha \trasp v_k A v_k + \trasp v_k (A x_k - b) = 0
\]
e quindi ponendo $r_k = b - Ax_k$ si ottiene
\[
 \alpha_k = \frac{\trasp v_k (-Ax_k + b)}{\trasp v_k A v_k} = \frac{\trasp v_k r_k}{\trasp v_k A v_k}
\]
Osserviamo in particolare che tenere traccia del valore di $r_k$ è utile per decidere quando
fermare il metodo. $||r_k|| = ||Ax_k - b||$ è indice di quanto siamo ``distanti'' dalla soluzione.
 
Si è verificato che questo metodo è convergente per ogni scelta di $x_0$ però ci si è accorti che
non è la scelta migliore della direzione di decrescita.
 
\subsection{Il metodo del gradiente coniugato}
Dopo l'introduzione del metodo del gradiente ottimo si è studiata un altro metodo di scelta, a cui
dobbiamo però premettere alcune definizioni.
 
\begin{de} Data una matrice $A \in \mat{\R}{n}$ una $n$-upla di vettori $(p_1 \ldots p_n)$ di $\R^n$
 si dicono \emph{$A$-coniugati} se
 \[
  \left[ \begin{array}{ccc}
          \quad \trasp p_1 \quad \\
          \quad \vdots \quad \\
          \quad \trasp p_n \quad \\
         \end{array} \right]
 A
 \left[ \begin{array}{ccc}
	\multirow{3}{*}{$p_1$} & \multirow{3}{*}{$\ldots$} & \multirow{3}{*}{$p_n$} \\
        & & \\
        & &
        \end{array} \right]
 = D
 \]
 dove $D$ è una matrice diagonale.
\end{de}
 
\begin{os} \label{os:coniugli}
 Se abbiamo una $n$-upla di vettore $A$-coniugati ed $A$ è invertibile
 possiamo facilmente dire che sono linearmente
 indipendenti.
\end{os}
 
Vorremmo ora impostare l'iterazione in modo che i vettori $v_0, v_1, \ldots$ che definiscono le
direzioni di descrescita siano $A$-coniugati. Osserviamo in particolare che la condizione di
lineare indipendenza ci permette di dire che non possiamo fare più di $n$ iterazioni dove $n$ è
la dimensione della matrice $A$. Scopriremo, per fortuna, che questo è solo apparentemente un
problema.
 
Dobbiamo ora trovare un metodo dati $v_1, \ldots, v_{k-1}$ vettori $A$-coniugati per determinare
$v_k$ tale che sia $A$-coniugato con tutti gli altri. Poniamo $v_k = r_k + \beta_k v_{k-1}$;
per avere la condizione di ortogonalità\footnote{secondo il prodotto scalare definito da $A$}
dovremo avere
\[
 (\trasp r_k + \beta_k\trasp v_{k-1}) A v_{k-1} = 0
\]
e quindi
\[
 \beta_k = \frac{\trasp r_k A v_{k-1}}{\trasp v_{k-1} A v_{k-1}}
\]
Si può verificare che questa condizione non è solo necessaria ma anche sufficiente. Un successione di vettori
scelti in questo modo è infatti forzatamente $A$-coniugata\footnote{questo risultato non verrà dimostrato qui.}
 
Avendo risolto il problema del trovare effettivamente la successione, ci chiediamo come affrontare
il limite al numero di iterazioni che ci è indicato Osservazione~\ref{os:coniugli}. Ci aiuta il
seguente
\begin{te}
 Siano $\{x_k\}_{k=1,\ldots,n}$ una successione di vettori che rispetti le condizioni osservate precedentemente con $x_0 = 0$;
 Allora per ogni $k$ si ha che
 \[
  \Phi(x_k) = \min_{x \in \Span{v_0, \ldots, v_k-1}}(\Phi(x))
 \]
 ed in particolare dopo $n$ iterazioni $\Phi(x_n)$ è il minimo assoluto di $\Phi(x)$.
\end{te}
Questo teorema ci dice, in particolare, che questo metodo iterativo è in realtà un metodo diretto, ovvero
è convergente in un numero finito di passi per ogni scelta di punto iniziale.
 
Poniamo ora $e_k = x_k - x$ l'errore al passo $k-esimo$. Il teorema sopra ci dice che $e_n = 0$ ma ci chiediamo
come si comporta $e_k$ mentre $k$ si avvicina $n$. Non dobbiamo infatti dimenticare che questo è un metodo
iterativo e quindi in molti casi pratici saremo interessati a fermarci prima di avere la convergenza totale.
Esiste un altro teorema che ci dà un risultato piuttosto preciso sulla velocità di convergenza:
\begin{te}
 Sia $||\cdot||_A$ la norma indotta dalla matrice $A$\footnote{ricordiamo che questa è definita nel seguente modo
 $||x||_A := \sqrt{\trasp x A x}$}
 ; allora per ogni $k$ e per ogni scelta iniziale di $x_0$ si
 ha
 % TODO: Questa maggiorazione è sbagliata. Il Demmel ne propone un'altra, provare a vedere se riusciamo
 % a portarla a qualcosa di umano.
 \[
  ||e_k||_A \leq \left( \frac{2\sqrt{\cond{A}}}{\sqrt{\cond{A}} - 1} \right)^{k} ||e_0||_A
 \]
\end{te}
Ancora una volta la velocità di convergenza dipende dal condizionamento della matrice del problema. Per
matrici con un piccolo condizionamento questa è molto veloce, mentre per matrici con condizionamento
più grande potrebbe diventare lenta.
 
\subsection{Precondizionamento}
Supponiamo di avere una matrice $A$ definita positiva che individua un problema mal condizionato, in cui la velocità
di convergenza del metodo del gradiente coniugato è molto lenta.
 
Sappiamo dai risultati di Analisi Numerica che non possiamo riuscire a risolvere con precisione
un sistema mal condizionato; ci chiediamo però se sia possibile perlomeno risolverlo velocemente,
pur con la consapevolezza che i risultati saranno molto imprecisi.
 
La risposta è sì e l'idea è usare un \emph{precondizionamento}, ovvero analizzare un altro problema
che si ottiene dal primo moltiplicandolo a destra e/o a sinistra per delle altre matrici rendendolo
ben condizionato. Ovviamente i risultati finali risentiranno del cattivo condizionamento del problema
iniziale, ma la velocità di convergenza sarà elevata.
 
Dato il sistema $Ax = b$ consideriamo il seguente
\[
 LA\trasp L {\trasp L}^{-1} x = L b
\]
che è equivalente. Osserviamo poi che $LA\trasp L$ è simile a $L^{-1}( L A \trasp L) L$ e quindi
a $A\trasp L L$. Ricordiamo ora che il condizionamento della matrice è dato da $\frac{\lambda_{max}}{\lambda_{min}}$
e cerchiamo quindi una matrice $M = \trasp L L $ tale che gli autovalori di $A \trasp LL $ siano molto vicini
\footnote{e quindi il loro rapporto sia molto vicino ad $1$, che è il valore minimo in cui possiamo sperare
per il condizionamento di una matrice}.
 
Osserviamo che se fosse possibile, ad esempio, scegliere $M = A^{-1}$ allora avremmo la situazione migliore possibile.
Evidentemente però questa non è un opzione, perché se fossimo a conoscenza di $A^{-1}$ avremmo già
completamente risolto il nostro problema. In ogni caso, però, un buon precondizionamento si ottiene cercando
di approssimare un'inversa di $A$.
 
Non ci occuperemo in questo corso di tecniche di precondizionamento che sono varie e a volte complesse.
Sottolineamo solo alcune problematiche che potrebbero nascere usando questo approccio
\begin{enumerate}[(a)]
 \item Per come abbiamo definito le matrici in gioco, $M$ dovrebbe essere definita positiva; in realtà
 esiste un modo per ovviare a questo problema, ma non verrà esposto qui;
 \item Una volta trovata $M$ (che potrebbe essere un'operazione complicata) dovremmo trovare anche
 $L$ e quindi fare una fattorizzazione; ancora una volta, esiste un modo per ignorare il fatto che esista
 la matrice $L$ ed usare solo $M$;
 \item Se $A$ è una matrice strutturata, saremo probabilmente interessati a mantenerne la struttura. Questo si può
 fare accontentandosi di approssimazioni dell'inversa piuttosto vaghe, ma poco invasive (come ad esempio una matrice diagonale);
\end{enumerate}
 
\section{Matrici strutturate}
Uno dei punti cardine nell'applicazione dei metodi per la soluzione dei sistemi lineari sarà riuscire a
sfruttare la struttura delle matrici con cui dovremo operare.
 
\subsection{Le matrici di Toeplitz} \label{subsec:toeplitz}
Molte delle matrici con cui avremo a che fare risolvendo sistemi lineari in pratica sono matrici di Toeplitz.
\begin{de}
 Sia $A$ una matrice $n \times n$; $A$ si dice di \emph{Toeplitz} se le sue diagonali sono costanti, ovvero se per
 ogni $i,j$ e per ogni $k \in \Z$ per cui ha senso
 \[
  a_{ij} = a_{i+k,j+k}
 \]
\end{de}
Queste matrici hanno una struttura molto particolare, ed esiste un modo piuttosto comodo di effettuare il
prodotto matrice per vettore. Consideriamo il caso seguente con una matrice di Toeplitz triangolare inferiore
\[
 \left[ \begin{array}{cccc}
         t_0 &  &  &  \\
         t_1 & \ddots &  &  \\
          \vdots & \ddots & \ddots &  \\
          t_n & \cdots & t_1 & t_0 \\
        \end{array} \right]
 \left[ \begin{array}{c}
         p_0 \\
         p_1 \\
         \vdots  \\
         p_n
        \end{array} \right]
= \left[ \begin{array}{l}
          t_0p_0 \\
          t_1p_0 + t_0p_1 \\
          \vdots\\
          t_np_0 + t_{n-1}p_1 + \ldots + t_0p_n
         \end{array} \right]
\]
Si osserva che il vettore che si ottiene ha i coefficienti che sono quelli del prodotto di questi due polinomi
\[
 \left\{ \begin{array}{lll}
          t(x) &=& t_0 + t_1 z + \ldots + t_n z^n \\
          p(x) &=& p_0 + p_1 z + \ldots + p_n z^n
         \end{array} \right.
\]
Possiamo quindi calcolare il prodotto matrice vettore nello stesso modo in cui calcoleremmo i coefficienti del
polinomio prodotto, ovvero con la trasformata discreta di Fourier.
 
Avremo un costo delle moltiplicazione $O(n\log(n))$\footnote{utilizzando la \fft.} e quindi
un costo complessivo del metodo del gradiente coniugato di $O(n^2\log(n))$.
 
% TODO: Inserire gli esempi delle applicazioni del metodo del gradiente a qualche caso particolare
% di matrici, come ad esempio le matrici elementari e le matrici con nugoli di autovalori appiccicati.
 
 
\subsection{Matrici di Toeplitz tridiagonali simmetriche}
Vorremmo ora mostrare un'analisi di un caso paricolare, ovvero dei sistemi lineari con
una matrice di Toeplitz simmetrica tridiagonale. \\
Questo è ad esempio il caso che discretizza il problema differenziale $\triangle u = f$ nel caso di
$u$ in una variabile reale. Le conclusione che otteremo su questo caso particolare ci permetteranno poi
di analizzare anche i casi in più variabili.
 
Supponiamo ora di avere una qualsiasi matrice $T$ tridiagonale simmetrica di Toeplitz di dimensione $n$; chiamiamo $a$ gli
elementi sulla diagonale e $b$ quelli sulla sotto e sopradiagonale\footnote{in effetti queste due scelte
individuano completamente la matrice di cui stiamo parlando.}. \\
Siamo interessati a studiare le proprietà spettrali di questa matrice per poter dare una stima
del suo condizionamento in norma 2, che come abbiamo visto influenza la convergenza del metodo del
gradiente coniugato.
 
Osserviamo ceh se $\lambda$ è un autovalore per $T$ allora la matrice $T - \lambda I$ deve essere singolare
e in particolare deve esistere una soluzione non banale del sistema lineare
\[
 (T - \lambda I)x = 0  \iff  \left[ \begin{array}{ccccc}
         a - \lambda & b &  & & \\
         b & a -\lambda & b &  & \\
         & \ddots & \ddots & \ddots & \\
         & & b & a - \lambda & b \\
         & & & b & a
        \end{array} \right]
 \left[ \begin{array}{c}
         x_1 \\
         x_2 \\
         \vdots \\
         x_{n-1} \\
         x_n
        \end{array} \right] =
\left[ \begin{array}{c}
        0 \\
        0 \\
        \vdots \\
        0 \\
        0
       \end{array} \right]
\]
Preso un qualsiasi $j$ compreso fra $2$ e $n-1$ possiamo scrivere la relazione sopra come
\[
 bx_{j-1} + (a- \lambda)x_j + bx_{j+1} = 0
\]
Ponendo $x_0 = x_{n+1} = 0$ la relazione vale anche per $j = 1$ e $j=n$. Apparentemente non abbiamo
chiarito molto e trovare una soluzione esplicita sembra complicato. Possiamo però provare a porre
$x_j = x^{j}$, e vedere se riusciamo a trovare una soluzione particolare di questo tipo. Sostituendo
nelle relazioni sopra (per $j$ fra $2$ e $n-1$) si ottiene
\[
 bx^{j-1} + (a- \lambda)x^j + bx^{j+1} = 0
\]
e ricordando che l'autovettore non può essere nullo e quindi $x^j \neq 0$ per ogni $j$, questa è soddisfatta
se e solo se
\[
 b + (a-\lambda)x + bx^2 = 0 \iff 1 + \frac{a - \lambda}{b} + x^2 = 0
\]
Passiamo ora ad analizzare il caso che ci interessa, ovvero $a = 2$ e $b = 1$\footnote{ovvero la matrice che discretizza
il problema differenziale di Laplace.}; per il terzo teorema di Gerschgorin sappiamo che $|\lambda - 2| < 2$ e
quindi
\[
 \left| \frac{a - \lambda}{b} \right|< 2
\]
Possiamo allora porre $\frac{a - \lambda}{b} = -2 \cos\theta$ per $\theta \in (0,\pi)$ e otteniamo l'equazione
\[
 1 - 2x\cdot\cos \theta + x^2 = 0
\]
Risolvendola otteniamo
\[
  x_{1,2} = \cos \theta \pm \sqrt{\cos^2\theta - 1} = \cos \theta \pm i \cdot \sin \theta = e^{\pm i \theta}
\]
Abbiamo quindi ottenuto due soluzioni che andrebbe bene per tutte le $j = 2, \ldots, n-1$ ma nessuna delle due
soddisfa le condizione per $j = 0$ e $j = n$. Osserviamo però che una qualsiasi combinazione lineare
$\alpha x_1 + \beta x_2$ soddisfa le condizione interne, e cerchiamo quindi di determinare $\alpha$ e $\beta$
in modo che anche le condizioni al contorno siano soddisfatte. Si ottiene
\[
 x_0 = 0 = \alpha + \beta
\]
e quindi $\beta = -\alpha$; ponendo poi $j = n+1$ si ottiene
\[
 x_{n+1} = 0 = \alpha e^{i\cdot(n+1) \theta} - \alpha e^{-i \cdot (n+1) \theta} = \alpha( 2 i \sin ((n+1)\theta) )
\]
e quindi $\theta_k = \frac{k\pi}{n+1}$ con $k = 1, \ldots, n$. Abbiamo trovato quindi $n$ autovettori distinti
e siamo in grado di calcolare i relativi autovalori ricordando che
\[
 \frac{a - \lambda_k}{b} = -2 \cos \theta_k = -2 \cos (\frac{k\pi}{n+1}) \Rightarrow \lambda_k = a + 2b\cos \theta_k
\]
Se costruiamo la matrice degli autovettori
\[
 U = \left[ \begin{array}{c|c|c|c}
             \multirow{3}{*}{$x_1$} & \multirow{3}{*}{$x_2$} & \multirow{3}{*}{$\ldots$} & \multirow{3}{*}{$x_n$} \\
             & & & \\
             & & &
            \end{array} \right]
\]
possiamo osservare che $\trasp UU = D$ con $D = \gamma I$. Abbiamo quindi che $U$ è quasi unitaria, in particolare
$\frac{1}{\gamma} U$ lo è. Inoltre possiamo osservare che gli elementi di $U$ non dipendono da $a$ e da $b$ e che
$u_{ij} = \sin(\frac{ij\pi}{n+1}) = u_{ji}$ e quindi $U$ è simmetrica. In altre parole abbiamo una decomposizione
spettrale $T = UDU$ dove tutta l'informazione sulla matrice è contenuta nella parte diagonale.
 
Osserviamo infine che $D = \diag{a + 2b\cos(\theta_1) , \ldots, a + 2b \cos(\theta_n)}$ e quindi l'autovalore più
piccolo è $a + 2b\cos(\theta_1)$.