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|\ge |{\lambda }_3|\ge \dots \ge |{\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 ${𝒙}_1,\dots ,{𝒙}_n$ autovettori linearmente indipendenti, che formano una base.

Consideriamo un vettore

${𝒛}_0=\sum _{j=1}^n{\alpha }_j{𝒙}_j,\text{rappresentazione }\circ$

rappresentato nella base degli autovettori. Richiedo che nella scrittura di ${𝒛}_0$ ci sia un contributo dell'autovettore ${𝒙}_1$ relativo all'autovalore ${\lambda }_1$, cioè ${\alpha }_1\ne 0$.

Genero una successione di vettori ${𝒛}_1,\dots ,{𝒛}_k$ tali che

${𝒛}_k=A\ast {𝒛}_{k-1}=A^k\ast {𝒛}_0\text{ da cui metodo della potenza}$

Allora sostituendo ${𝒛}_0$ la rappresentazione $\circ$ ottengo:

${𝒛}_k=A^k[\sum _{j=1}^n{\alpha }_j{𝒙}_j]$

e siccome gli ${𝒙}_j$ sono autovettori possiamo scrivere:

${𝒛}_k=\sum _{j=1}^n{\lambda }_j^k{\alpha }_j{𝒙}_j$

${𝒛}_k={\lambda }_1^k{\alpha }_1{𝒙}_1+\sum _{j=2}^n{\lambda }_j^k{\alpha }_j{𝒙}_j$

e raccogliendo ${\lambda }_1^k$:

${𝒛}_k={\lambda }_1^k[{\alpha }_1{𝒙}_1+\sum _{j=2}^n(\frac{{\lambda }_j}{{\lambda }_1}{)}^k{\alpha }_j{𝒙}_j],\text{rappresentazione }\mathrm{♡}$

${\lambda }_j<{\lambda }_1⟶(\frac{{\lambda }_j}{{\lambda }_1}{)}^k\to 0⟺k\to \mathrm{\infty }$, allora per $k\to \mathrm{\infty }$, ${𝒛}_k\parallel {𝒙}_1$.

Definisco il quoziente di Rayleigh:

${\sigma }_k=\frac{{𝒛}_k^h\ast A\ast {𝒛}_k}{{𝒛}_k^h{𝒛}_k}$

Verifichiamo che ${\sigma }_k\to {\lambda }_1$. Sostituisco l'espressione $\mathrm{♡}$ di ${𝒛}_k$ nel ${\sigma }_k$ che abbiamo definito:

${\sigma }_k=\frac{[{\overline{\lambda }}_1^k({\overline{\alpha }}_1{𝒙}_1^h+\sum _{j=2}^n{\overline{(\frac{{\lambda }_j}{{\lambda }_1})}}^k{\overline{\alpha }}_j{𝒙}_j^h]\ast A\ast [{\lambda }_1^k({𝒙}_1{\alpha }_1+\sum _{j=2}^n(\frac{{\lambda }_j}{{\lambda }_1}{)}^k{\alpha }_j{𝒙}_j]}{[{\overline{\lambda }}_1^k({\overline{\alpha }}_1{𝒙}_1^h+\sum _{j=2}^n{\overline{(\frac{{\lambda }_j}{{\lambda }_1})}}^k{\overline{\alpha }}_j{𝒙}_j^h]\ast [{\lambda }_1^k({\alpha }_1{𝒙}_1+\sum _{j=2}^n(\frac{{\lambda }_j}{{\lambda }_1}{)}^k{\alpha }_j{𝒙}_j]}$

Per $k\to \mathrm{\infty }$ i termini $(\frac{{\lambda }_j}{{\lambda }_1}{)}^k\to 0$, mentre i termini ${\overline{\lambda }}_1^k\ast {\lambda }_1^k$ si semplificano e rimane:

$\underset{k\to \mathrm{\infty }}{lim}{\sigma }_k=$

$=\frac{{\overline{\alpha }}_1{\alpha }_1{𝒙}_1^{h}A{𝒙}_1}{{\overline{\alpha }}_1{\alpha }_1{𝒙}_1^h{𝒙}_1}$

$=\frac{{𝒙}_1^h{\lambda }_1{𝒙}_1}{{𝒙}_1^h{𝒙}_1}={\lambda }_1$

Sul calcolatore il problema compare nella divisione per ${𝒙}_k^h{𝒙}_k=|{𝒙}_k|$, perché si rischia di andare nella barriera di overflow o di underflow. Per risolvere il problema, invece di considerare ${𝒛}_k$, lo si normalizza in norma 2 cioè si considera:

${𝒚}_k=\frac{{𝒛}_k}{|{𝒛}_k|}$

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

Implementazione dell'algoritmo

Dato ${𝒛}_0$, eseguo le seguenti operazioni:

${𝒚}_k=\frac{{𝒛}_k}{|{𝒛}_k|}$

${𝒛}_{k+1}=A{𝒚}_k$

$\sigma ={𝒚}_k^h{𝒛}_{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\ast (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\ne 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}|<\epsilon$

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

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

Proposizione 3.1

Supponiamo di avere $A\in {ℂ}_{n,n}$ normale, consideriamo $𝒙\in {ℂ}^n\setminus 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 )|\le \frac{|f(A)𝒙{|}_2}{|𝒙{|}_2}$

Test 1: Impongo la condizione

$|{𝒛}_{k+1}-{\sigma }_k{𝒚}_k|\le \epsilon$

che permette di avere un controllo sull'errore assoluto.

Infatti, se scegliamo la funzione $f(\lambda )=\lambda -{\sigma }_k$ (errore assoluto) e $𝒙={𝒚}_k$, allora, per la proposizione:

$|f(\lambda )|=|\lambda -{\sigma }_k|\le \frac{|(A-{\sigma }_{k}I){𝒚}_k|}{|{𝒚}_k{|}_2}$

$|{𝒚}_k|=1$ per costruzione, allora

$|\lambda -{\sigma }_k|\le |A\ast {𝒚}_k-{\sigma }_{k}I{𝒚}_k|=|{𝒛}_{k+1}-{\sigma }_k{𝒚}_k|$

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

$\frac{|{𝒛}_{k+1}-{\sigma }_k{𝒚}_k{|}_2}{|{𝒛}_{k+1}|}\le \epsilon$

permette di controllare l'errore relativo.

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

$\frac{|\lambda -{\sigma }_k|}{\lambda }\le \frac{|(A-{\sigma }_{k}I)A^{-1}{𝒛}_{k+1}{|}_2}{|{𝒛}_{k+1}|}$

$=\frac{|{𝒛}_{k+1}-{\sigma }_k{A}^{-1}{𝒛}_{k+1}|}{|𝒛k+1{|}_2}$

$=\frac{|{𝒛}_{k+1}-{\sigma }_k{𝒚}_k|}{|{𝒛}_{k+1}{|}_2}$