Calcolo numerico/Metodo della potenza

Autovalore di modulo massimo

Supponiamo che l'autovalore di modulo massimo $|\lambda _{1}|$ sia strettamente maggiore dell'autovalore successivo, cioè

$|\lambda _{1}|>|\lambda _{2}|\geq |\lambda _{3}|\geq \dots \geq |\lambda _{n}|$

si esclude quindi il caso in cui $\lambda _{1}$ è un autovalore complesso (e avrebbe modulo uguale a quello del suo coniugato) o in cui $\lambda _{1}$ ha molteplicità maggiore di 1 (non sarebbe strettamente maggiore di se stesso).

Supponiamo che esistano ${\boldsymbol {x}}_{1},\dots ,{\boldsymbol {x}}_{n}$ autovettori linearmente indipendenti, che formano una base.

Consideriamo un vettore

${\boldsymbol {z}}_{0}=\sum _{j=1}^{n}\alpha _{j}{\boldsymbol {x}}_{j},\quad {\hbox{rappresentazione }}\circ$

rappresentato nella base degli autovettori. Richiedo che nella scrittura di ${\boldsymbol {z}}_{0}$ ci sia un contributo dell'autovettore ${\boldsymbol {x}}_{1}$ relativo all'autovalore $\lambda _{1}$, cioè $\alpha _{1}\neq 0$.

Genero una successione di vettori ${\boldsymbol {z}}_{1},\dots ,{\boldsymbol {z}}_{k}$ tali che

${\boldsymbol {z}}_{k}=A*{\boldsymbol {z}}_{k-1}=A^{k}*{\boldsymbol {z}}_{0}\quad {\hbox{ da cui metodo della potenza}}$

Allora sostituendo ${\boldsymbol {z}}_{0}$ la rappresentazione $\circ$ ottengo:

${\boldsymbol {z}}_{k}=A^{k}[\sum _{j=1}^{n}\alpha _{j}{\boldsymbol {x}}_{j}]$

e siccome gli ${\boldsymbol {x}}_{j}$ sono autovettori possiamo scrivere:

${\boldsymbol {z}}_{k}=\sum _{j=1}^{n}\lambda _{j}^{k}\alpha _{j}{\boldsymbol {x}}_{j}$

${\boldsymbol {z}}_{k}=\lambda _{1}^{k}\alpha _{1}{\boldsymbol {x}}_{1}+\sum _{j=2}^{n}\lambda _{j}^{k}\alpha _{j}{\boldsymbol {x}}_{j}$

e raccogliendo $\lambda _{1}^{k}$:

${\boldsymbol {z}}_{k}=\lambda _{1}^{k}[\alpha _{1}{\boldsymbol {x}}_{1}+\sum _{j=2}^{n}({\frac {\lambda _{j}}{\lambda _{1}}})^{k}\alpha _{j}{\boldsymbol {x}}_{j}],\quad {\hbox{rappresentazione }}\heartsuit$

$\lambda _{j}<\lambda _{1}\,\longrightarrow \,({\frac {\lambda _{j}}{\lambda _{1}}})^{k}\to 0\iff k\to \infty$, allora per $k\to \infty$, ${\boldsymbol {z}}_{k}\parallel {\boldsymbol {x}}_{1}$.

Definisco il quoziente di Rayleigh:

$\sigma _{k}={\frac {{\boldsymbol {z}}_{k}^{h}*A*{\boldsymbol {z}}_{k}}{{\boldsymbol {z}}_{k}^{h}{\boldsymbol {z}}_{k}}}$

Verifichiamo che $\sigma _{k}\to \lambda _{1}$. Sostituisco l'espressione $\heartsuit$ di ${\boldsymbol {z}}_{k}$ nel $\sigma _{k}$ che abbiamo definito:

$\sigma _{k}={\frac {[{\bar {\lambda }}_{1}^{k}({\bar {\alpha }}_{1}{\boldsymbol {x}}_{1}^{h}+\sum _{j=2}^{n}{\bar {({\frac {\lambda _{j}}{\lambda _{1}}})}}^{k}{\bar {\alpha }}_{j}{\boldsymbol {x}}_{j}^{h}]*A*[\lambda _{1}^{k}({\boldsymbol {x}}_{1}\alpha _{1}+\sum _{j=2}^{n}({\frac {\lambda _{j}}{\lambda _{1}}})^{k}\alpha _{j}{\boldsymbol {x}}_{j}]}{[{\bar {\lambda }}_{1}^{k}({\bar {\alpha }}_{1}{\boldsymbol {x}}_{1}^{h}+\sum _{j=2}^{n}{\bar {({\frac {\lambda _{j}}{\lambda _{1}}})}}^{k}{\bar {\alpha }}_{j}{\boldsymbol {x}}_{j}^{h}]*[\lambda _{1}^{k}(\alpha _{1}{\boldsymbol {x}}_{1}+\sum _{j=2}^{n}({\frac {\lambda _{j}}{\lambda _{1}}})^{k}\alpha _{j}{\boldsymbol {x}}_{j}]}}$

Per $k\to \infty$ i termini $({\frac {\lambda _{j}}{\lambda _{1}}})^{k}\to 0$, mentre i termini ${\bar {\lambda }}_{1}^{k}*\lambda _{1}^{k}$ si semplificano e rimane:

$\lim _{k\to \infty }\sigma _{k}=$

$={\frac {{\bar {\alpha }}_{1}\alpha _{1}{\boldsymbol {x}}_{1}^{h}A{\boldsymbol {x}}_{1}}{{\bar {\alpha }}_{1}\alpha _{1}{\boldsymbol {x}}_{1}^{h}{\boldsymbol {x}}_{1}}}$

$={\frac {{\boldsymbol {x}}_{1}^{h}\lambda _{1}{\boldsymbol {x}}_{1}}{{\boldsymbol {x}}_{1}^{h}{\boldsymbol {x}}_{1}}}=\lambda _{1}$

Sul calcolatore il problema compare nella divisione per ${\boldsymbol {x}}_{k}^{h}{\boldsymbol {x}}_{k}=\vert {\boldsymbol {x}}_{k}\vert$, perché si rischia di andare nella barriera di overflow o di underflow. Per risolvere il problema, invece di considerare ${\boldsymbol {z}}_{k}$, lo si normalizza in norma 2 cioè si considera:

${\boldsymbol {y}}_{k}={\frac {{\boldsymbol {z}}_{k}}{\vert {\boldsymbol {z}}_{k}\vert }}$

e in questo modo nel quoziente di Rayleigh si divide per un vettore di norma 1.

Implementazione dell'algoritmo

Dato ${\boldsymbol {z}}_{0}$, eseguo le seguenti operazioni:

${\boldsymbol {y}}_{k}={\frac {{\boldsymbol {z}}_{k}}{\vert {\boldsymbol {z}}_{k}\vert }}$

${\boldsymbol {z}}_{k+1}=A{\boldsymbol {y}}_{k}$

$\sigma ={\boldsymbol {y}}_{k}^{h}{\boldsymbol {z}}_{k+1}$

e se una certa condizione imposta nel test d'arresto non è verificata, si procede.

Il costo è dato da $n^{2}$ operazioni del prodotto matrice-vettore, e $3n$ operazioni per i prodotti scalari nelle norme.

$\sigma _{k}=\lambda _{1}*(1+O({\frac {\lambda _{2}}{\lambda _{1}}})^{k}))$, cioè la velocità con cui $\sigma _{k}$ converge a $\lambda _{1}$ è come il rapporto tra ${\frac {\lambda _{2}}{\lambda _{1}}}$, e se i due autovalori sono ben separati il metodo converge rapidamente.

Il metodo converge anche se l'autovalore di modulo massimo non ha molteplicità algebrica uguale a 1, in particolare:

1. $\lambda _{1}=\lambda _{2}$. I due (o $r$) autovalori dominanti sono coincidenti. In questo caso il metodo risulta ancora convergente, ma $x_{k}$ converge a un vettore nel sottospazio generato da $x_{1}$e $x_{2}$;
2. $\lambda _{1}=-\lambda _{2}$. I due autovalori dominanti sono opposti. In questo caso l'autovalore di modulo massimo può essere approssimato applicando il metodo delle potenze alla matrice $A^{2}$, in modo che $\lambda _{1}^{2}=\lambda _{2}^{2}$ e l'analisi ricada nel punto precedente;
3. i due autovalori dominanti sono complessi coniugati. In questo caso si hanno delle oscillazioni non smorzate nella sequenza dei vettori $x_{k}$ e il metodo delle potenze non converge.

$\alpha _{1}\neq 0$ non è un'ipotesi essenziale se si lavora in aritmetica floating point.

Test di arresto

Test dell'incremento: richiediamo che

$|\sigma _{k}-\sigma _{k-1}|<\varepsilon$

e il metodo prosegue fino a quando questa condizione non è verificata.

Nel caso in cui $A$ è una matrice normale, cioè $AA_{h}=A_{h}A$, possiamo considerare altri due test d'arresto.

Proposizione 3.1

Supponiamo di avere $A\in \mathbb {C} _{n,n}$ normale, consideriamo ${\boldsymbol {x}}\in \mathbb {C} ^{n}\smallsetminus 0$, e $f$ razionale definita su un sottoinsieme del piano complesso che contiene tutti gli autovalori di $A$. Allora esiste $\lambda \in A$ tale che

$|f(\lambda )|\leq {\frac {\vert f(A){\boldsymbol {x}}\vert _{2}}{\vert {\boldsymbol {x}}\vert _{2}}}$

Test 1: Impongo la condizione

$\vert {\boldsymbol {z}}_{k+1}-\sigma _{k}{\boldsymbol {y}}_{k}\vert \leq \varepsilon$

che permette di avere un controllo sull'errore assoluto.

Infatti, se scegliamo la funzione $f(\lambda )=\lambda -\sigma _{k}$ (errore assoluto) e ${\boldsymbol {x}}={\boldsymbol {y}}_{k}$, allora, per la proposizione:

$|f(\lambda )|=|\lambda -\sigma _{k}|\leq {\frac {\vert (A-\sigma _{k}I){\boldsymbol {y}}_{k}\vert }{\vert {\boldsymbol {y}}_{k}\vert _{2}}}$

$\vert {\boldsymbol {y}}_{k}\vert =1$ per costruzione, allora

$|\lambda -\sigma _{k}|\leq |A*{\boldsymbol {y}}_{k}-\sigma _{k}I{\boldsymbol {y}}_{k}|=\vert {\boldsymbol {z}}_{k+1}-\sigma _{k}{\boldsymbol {y}}_{k}\vert$

Test 2: se $A$ è normale e non singolare, allora la condizione

${\frac {\vert {\boldsymbol {z}}_{k+1}-\sigma _{k}{\boldsymbol {y}}_{k}\vert _{2}}{\vert {\boldsymbol {z}}_{k+1}\vert }}\leq \varepsilon$

permette di controllare l'errore relativo.

Scegliamo $f(\lambda )={\frac {\lambda -\sigma _{k}}{\lambda }}$ (errore relativo), e poniamo ${\boldsymbol {x}}={\boldsymbol {z}}_{k+1}$, applichiamo poi la proposizione:

${\frac {|\lambda -\sigma _{k}|}{\lambda }}\leq {\frac {\vert (A-\sigma _{k}I)A^{-1}{\boldsymbol {z}}_{k+1}\vert _{2}}{\vert {\boldsymbol {z}}_{k+1}\vert }}$

$={\frac {\vert {\boldsymbol {z}}_{k+1}-\sigma _{k}A^{-1}{\boldsymbol {z}}_{k+1}\vert }{\vert {\boldsymbol {z}}{k+1}\vert _{2}}}$

$={\frac {\vert {\boldsymbol {z}}_{k+1}-\sigma _{k}{\boldsymbol {y}}_{k}\vert }{\vert {\boldsymbol {z}}_{k+1}\vert _{2}}}$