From 2360e79e9eef0ae866ddcf0e17a62f60d974db1a Mon Sep 17 00:00:00 2001 From: Leonardo Robol Date: Fri, 23 Oct 2009 15:49:04 +0200 Subject: [PATCH] Lezione di Gemignani del 23 ottobre e modifiche nell'introduzione --- capitolo2.tex | 124 ++++++++++++++++++++++++++++++++++++++++++++++++++++++- introduzione.tex | 4 ++ 2 files changed, 127 insertions(+), 1 deletion(-) diff --git a/capitolo2.tex b/capitolo2.tex index 23c05e2..d31a7a3 100644 --- a/capitolo2.tex +++ b/capitolo2.tex @@ -531,4 +531,126 @@ Come ultima nota, ricordiamo che questo algoritmo è applicabile solo alle matri 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. \ No newline at end of file +robusta del metodo Divide et Impera sulle matrici di Hessenberg. + +\section{Il metodo delle potenze} +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 ottnere 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 raffinamente 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} + diff --git a/introduzione.tex b/introduzione.tex index 8df7a72..d74cd15 100644 --- a/introduzione.tex +++ b/introduzione.tex @@ -27,6 +27,10 @@ Questi appunti si dividono in due parti Durante il corso del prof. Bini abbiamo usato il linguaggio Fortrano 95 e quindi gli esempi saranno esposti con quel linguaggio. Sto comunque valutando l'opportunità di esporre anche qualche algoritmo in qualche linguaggio diverso. + \item[Parte III -- Appendice] Tutte le cose che in realtà non fanno strettamente parte del corso ma + che potrebbero essere utili per comprendere meglio quello che si trova scritto in questi appunti. Ci + si trova ad esempio, una breve introduzione al linguaggio Fortran che è usato nella parte II per + implementare gli algoritmi. \end{description} \section*{Notazione} -- 2.1.4