diff --git a/AppVibrazioni.tex b/AppVibrazioni.tex index 1f867ce..c35f1dc 100644 --- a/AppVibrazioni.tex +++ b/AppVibrazioni.tex @@ -164,7 +164,7 @@ soluzione per un certo $t$ ha un costo computazionale molto alto. \`E quindi ina per la valutazione in tutti i punti che interessano a noi per ottenere un suono bem campionato (circa 44100 al secondo). \end{os} -\section{Sistemi continui} +\section{Sistemi continui} \label{sec:vibsistcont} \subsection{Generalizzazione del caso delle $N$ masse} Supponiamo di avere una corda tesa fra due supporti fissi, che viene, in un'istante $t = 0$ messa in una certa posizione con una data velocità iniziale. Siamo interessati a studiare il successivo moto @@ -300,4 +300,62 @@ il cambio di base che la diagonalizza è ortogonale e quindi ha condizionamento è data dalla norma du $v$. Vogliamo ora mostrare che $\frac{||\tau||}{||v||}$ è limitato superiormente da una costante, in modo -da provare la convergenza della nostra approssimazione. +da provare la convergenza della nostra approssimazione. + +Consideriamo la seguente serie di eguaglianze +\[ + \sqrt{h} || \tau ||_2 = ( \sum_{i=1}^{n} v^{(4)}(\xi_i)^2 + v{(4)}(\eta_i)^2)^{\frac 1 2} \sqrt{h} + = ( \sum_{i=1}^{n} h v^{(4)}(\xi_i)^2 + h v{(4)}(\eta_i)^2)^{\frac 1 2} +\] +Osserviamo che l'ultima espressione può essere intesa, per $h \to 0$, come un'approssimazione dell'integrale +di Riemann di $v^2(x)$, ed essendo l'integrale di una funzione continua su un compatto è limitato +superiormente da una costante. Abbiamo quindi provato la tesi. + +Introduciamo ora un metodo iterativo che possiamo vedere come una estensione del metodo delle potenze +e che ci permetterà di studiare nel dettaglio questa situazione. + +\subsection{Il metodo delle iterazioni ortogonali di sottospazi} +Il metodo delle potenze come è stato analizzato nella Sezione~\ref{sec:metodopotenze} è basato sull'idea +di prendere un vettore e iterare su di lui l'applicazione individuata dalla matrice $A$ di cui si vogliono +calcolare gli autovalori. \\ +Osserviamo che questa idea è equivalente (anzi, forse è meglio rappresentata) all'iterare un sottospazio +vettoriale di dimensione $1$. Supponiamo ora di voler calcolare, invece che 1, più autovalori delle matrice, +precisamente $k$ (dove $k < n$ ed $n$ è la dimensione delle matrice) ognuno di molteplicità 1\footnote{% +Questo solo per semplicità di esposizione. Non è complicato generalizzare l'idea}. Una naturale estensione del metodo +sarebbe iterare un sottospazio vettoriale di dimensione $k$. Questo dovrebbe convergere ad un sottospazio di dimensione +$k$ uguale alla somma degli autospazi dei primi $k$ autovalori.\\ +Si può osservare subito un problema computazionale. Sebbene questo discorso sia teoricamente corretto, è facile prevedere +che il nostro spazio di dimensione $k$ convergerà verso l'autospazio relativo +all'autovalore dominante (che è di dimensione 1). Più precisamente, supponiamo di rappresentare lo spazio +con una matrice $x$ $n \times k$ le cui colonne siano una base dello spazio. Allora l'iterazione sarebbe +\[ +\left\{\begin{array}{lcl} + x_0 &=& x \\ + x_{k+1} &=& Ax_k +\end{array} \right. +\] +Per lo stesso ragionamento fatto nel caso del metodo delle potenze originale tutte le colonne di $A$ convergeranno +a multipli dello stesso vettore (quello dell'autovalore dominante). + +Un modo di risolvere questo problema potrebbe essere assicurarsi che durante l'iterazione tutte le colonne restino +ortogonali fra loro. Consideriamo la seguente iterazione +\[ + \left\{ \begin{array}{lcl} + x_0 &=& x \\ + T_{k} &=& A x_k \\ + Q_k R_k &=& T_k \: \text{(scomposizione $QR$ di $T_k$)} \\ + x_{k+1} &=& Q_k + \end{array} \right. +\] +Ricordando che la scomposizione QR di una matrice $n \times k$ con $k < n$ è tale che $Q$ è rettangolare +ed è $n \times k$ e $R$ è quadrata $k \times k$. In particolare si osserva che nella nostra iterazione +le colonne di $Q_k$ sono una +base ortonormale per lo spazio generato dalle colonne di $T_k$. +Osserviamo che calcolare esplicitamente $x_k$ non è realmente necessario e possiamo riformulare +l'iterazione. Poniamo +\[ + B_k = \trasp{Q_k} A Q_k \qquad \text{e} \qquad B_k = U_k D \trasp{U_k} +\] +dove $U_K D \trasp{U_k}$ è la decomposizione spettrale di $B_k$. Si osserva che $T_{k+1} = A Q_k = +Q_k B_k = $. +% TODO: Capire l'iterazione implicita. \ No newline at end of file diff --git a/CalcoloScientifico.tex b/CalcoloScientifico.tex index fc5ae85..8ab4b89 100644 --- a/CalcoloScientifico.tex +++ b/CalcoloScientifico.tex @@ -117,13 +117,17 @@ \newcommand{\deter}[1]{\mathrm{det}(#1)} \newcommand{\rk}[1]{\mathrm{rk}(#1)} \newcommand{\diag}[1]{\mathrm{diag}(#1)} +\newcommand{\Dim}[1]{\mathrm{dim}(#1)} %% %% I nomi di algoritmi e fattorizzazioni %% \newcommand{\svd}{\texttt{SVD}} \newcommand{\tsvd}{\texttt{TSVD}} -\newcommand{\Span}[1]{\mathrm{Span}(#1)} +\newcommand{\qr}{\texttt{QR}} +\newcommand{\matlab}{\texttt{Matlab}} +\newcommand{\lapack}{\texttt{Lapack}} +\newcommand{\Span}[1]{<(#1)>} \newcommand{\Ker}[1]{\mathrm{Ker}(#1)} \newcommand{\Exp}[1]{\mathrm{exp}{#1}} @@ -259,6 +263,9 @@ %% Decomposizione in valori singolari \include{capitolo3} +%% Risoluzione dei sistemi lineari +\include{capitolo4} + %% Comincia l'appendice, ovvero dove %% esponiamo gli algoritmi visti con Bini \part{Esercitazioni} diff --git a/capitolo2.tex b/capitolo2.tex index d31a7a3..dd628a1 100644 --- a/capitolo2.tex +++ b/capitolo2.tex @@ -533,7 +533,7 @@ 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} +\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 diff --git a/capitolo3.tex b/capitolo3.tex index 0a78ff7..1394092 100644 --- a/capitolo3.tex +++ b/capitolo3.tex @@ -276,3 +276,125 @@ il metodo di compressione, ma questo lo vedremo in seguito}. % TODO: Inserire un esempio di immagine compressa tramite SVD magari con varie opzioni per il rango +\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. + + + diff --git a/capitolo4.tex b/capitolo4.tex new file mode 100644 index 0000000..4151df3 --- /dev/null +++ b/capitolo4.tex @@ -0,0 +1,193 @@ +\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{Il metodo del gradiente coniugato} + +\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\footnote{nel senso di retta affine} $x_k + \Span{v_k}$: +\[ + 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 = v - Av_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 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} + Se ho $\{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 + \[ + ||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} +