From 6298d57eaab75ae46a1b61ae30532595971ec2c9 Mon Sep 17 00:00:00 2001 From: Leonardo Robol Date: Sat, 17 Oct 2009 17:11:07 +0200 Subject: [PATCH] Aggiunta lezione di gemignani del 16 ottobre. Ancora da controllare la parte sul Bulge Chasing di cui per altro pure io ho capito poco --- capitolo2.tex | 146 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 146 insertions(+) diff --git a/capitolo2.tex b/capitolo2.tex index 6236ad0..db11130 100644 --- a/capitolo2.tex +++ b/capitolo2.tex @@ -213,3 +213,149 @@ le considerazioni fatte sul primo passo valgono anche per i seguenti. In partico 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 quel caso 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 principalmente farebbe +perdere la simmetria\footnote{Per simmetria si intende 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 le cose continuano 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\protect\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.}} +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 Hessengerg. +Vogliamo quindi ora calcolare la fattorizzazione QR \textbf{mentre} calcoliamo $p_k(A_k)$ in modo da +rendere la risultante matrice in forma di Hessenberg e dimnuire 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. + +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. \ No newline at end of file -- 2.1.4