Calcolo di autovalori e autovettori

Alvise Sommariva

Università degli Studi di Padova

Dipartimento di Matematica

15 aprile 2014

Autovalori

Il problema del calcolo degli autovalori di una matrice quadrata A di ordine n consiste nel trovare gli n numeri (possibilmente complessi) λ tali che

Ax = λx, x ≠ 0  (1)

Si osservi che a seconda delle esigenze

▶ talvolta è richiesto solamente il calcolo di alcuni autovalori (ad esempio quelli di massimo modulo, per determinare lo spettro della matrice),

▶ talvolta si vogliono determinare tutti gli n autovalori in ℂ.

Autovalori

Per semplicità , dopo i teoremi di localizzazione di Gershgorin, mostreremo solo due metodi classici, uno per ognuna di queste classi, quello delle potenze e il metodo QR, rimandando per completezza alla monografia di Saad o a manuali di algebra lineare [2], [13].

Una interessante applicazione è l’algoritmo di PageRank [11], utilizzato da Google per fornire i risultati migliori tra i siti web relativamente a certe parole chiave ed in prima approssimazione basato sul calcolo di un autovettore relativo all’autovalore 1 (ad esempio via metodo delle potenze) di una matrice stocastica di dimensioni enormi.

Teoremi di Gershgorin

In questo paragrafo mostriamo tre teoremi di localizzazione di autovalori dovuti a Gershgorin (cf. [2, p.76], [15]).

Teorema

(Primo teorema di Gershgorin). Gli autovalori di una matrice A di ordine n sono tutti contenuti nell’unione dei cerchi di Gershgorin

K_{i} = {z ∈ ℂ : |z − a_{i,i}| ≤ ∑_{j=1,j≠i}^{n} |a_{i,j}|}

Vediamo quale esempio la matrice

A = ( \begin{matrix}15 & −2 & 2 \\ 1 & 10 & −3 \\ −2 & 1 & 0 \end{matrix} )  (2)

Il primo teorema di Gershgorin stabilisce che gli autovalori stanno nell’unione dei cerchi di Gershgorin

K_{1} = {z ∈ ℂ : |z − 15| ≤ | − 2| + |2| = 4}
K_{2} = {z ∈ ℂ : |z − 10| ≤ |1| + | − 3| = 4}
K_{3} = {z ∈ ℂ : |z − 0| ≤ | − 2| + |1| = 3}

Teorema

(Secondo teorema di Gershgorin). Se l’unione M_{1} di k cerchi di Gershgorin è disgiunta dall’unione M_{2} dei rimanenti n − k, allora k autovalori appartengono a M_{1} e n − k appartengono a M_{2}.

Applicando il secondo teorema di Gershgorin, dal confronto con la figura 6 abbiamo che 1 autovalore sta nel cerchio K_{3} mentre 2 stanno nell’unione dei cerchi K_{1}, K_{2}.

Alvise Sommariva  Calcolo di autovalori e autovettori  5/ 60

Teoremi di Gershgorin

Figura: Cerchi di Gershgorin della matrice A definita in (4)

Alvise Sommariva  Calcolo di autovalori e autovettori  6/ 60

Definizione

Una matrice di ordine n ≥ 2 è riducibile se esiste una matrice di permutazione Π e un intero k, 0 < k < n, tale che

B = ΠAΠ^{T} = ( \begin{matrix}A_{1,1} & A_{1,2} \\ 0 & A_{2,2} \end{matrix} )

in cui A_{1,1} ∈ ℂ^{k×k}, A_{2,2} ∈ ℂ^{(n−k)×(n−k)}.

Definizione

Una matrice si dice irriducibile se non è riducibile.

Nota.

Per verificare se una matrice A = (a_{i,j}) sia irriducibile (cf. [8]), ricordiamo che data una qualsiasi matrice, è possibile costruire un grafo avente come nodi gli indici della matrice. In particolare, il nodo i-esimo è connesso al nodo j-esimo se l’elemento a_{i,j} è diverso da 0.

Il grafo associato si dice fortemente connesso se per ogni coppia (i, j) posso raggiungere j a partire da i.

Una matrice è irriducibile se e solo se il grafo ad essa associata (detto di adiacenza) è fortemente connesso.

In altre parole, una matrice riducibile se e solo se il grafo di adiacenza ad esso associato non è fortemente connesso.

Teorema

(Terzo teorema di Gershgorin). Se la matrice di ordine n è irriducibile e un autovalore λ sta sulla frontiera dell’unione dei cerchi di Gershgorin, allora sta sulla frontiera di ogni cerchio di Gershgorin.

Esempio. La matrice

B = ( \begin{matrix}2 & −1 & 0 & 0 \\ −1 & 2 & −1 & 0 \\ 0 & −1 & 2 & −1 \\ 0 & 0 & −1 & 2 \end{matrix} )

è irriducibile, in quanto il grafico di adiacenza è fortemente connesso. Gli autovalori stanno nella disco centrato in (2, 0) e raggio 2 (primo teorema di Gershgorin), ma non sulla frontiera (terzo teorema di Gershgorin). Quindi è non singolare.

Alvise Sommariva  Calcolo di autovalori e autovettori  9/ 60

Vediamo ora in Matlab quali sono effettivamente gli autovalori. si ha

>> A=[15 −2 2 ; 1 10 −3; −2 1 0 ]

 A =

15  −2  2
1  10  −3
−2  1  0
>> e i g ( A )
ans =

0.5121 14.1026 10.3854

>>

a conferma di quanto stabilito dai primi due teoremi di Gershgorin.

Nota.

Ricordiamo che A è una matrice a coefficienti reali, allora A e A^{T} hanno gli stessi autovalori (cf. [2, p.47]) e quindi applicando i teoremi di Gershgorin alla matrice trasposta possiamo ottenere nuove localizzazioni degli autovalori.

Nel caso A sia a coefficienti complessi, se λ è un autovalore di A allora il suo coniugato \overline{λ} è autovalore della sua trasposta coniugata \overline{A}. Da qui si possono fare nuove stime degli autovalori di A.

Cosa possiamo dire relativamente agli autovalori di A se applichiamo i teoremi di Gershgorin ad A^{T} invece che ad A?

Metodo delle potenze

Il metodo delle potenze è stato suggerito nel 1913 da Muntz ed è particolarmente indicato per il calcolo dell’autovalore di massimo modulo di una matrice.

Sia A una matrice quadrata di ordine n con ▶ n autovettori x_{1}, …, x_{n} linearmente indipendenti, ▶ autovalori λ_{1}, …, λ_{n} tali che

|λ_{1}| > |λ_{2}| ≥ … ≥ |λ_{n}|.  (3)

A tal proposito ricordiamo (cf. [3], p. 951) i seguenti risultati.

▶ Una matrice A è diagonalizzabile se e solo se possiede n

autovettori linearmente indipendenti. ▶ Se tutti gli autovalori di A sono distinti la matrice è

diagonalizzabile; l’opposto è ovviamente falso (si pensi alla matrice identica).

▶ Una matrice simmetrica (hermitiana) è diagonalizzabile. L’opposto è ovviamente falso, visto che la matrice

A = ( \begin{matrix}15 & 0 \\ 1 & 10 \end{matrix} )  (4)

è diagonalizzabile visto che ha tutti gli autovalori distinti ma non è simmetrica.

Metodo delle potenze

Il metodo delle potenze è definito come segue. Sia t_{0} ∈ R^{n} definito da

t_{0} = ∑_{i=1}^{n}α_{i} x_{i}, α_{1} ≠ 0,

e si generi la successione

y_{0} = t_{0}  (5)
y_{k} = Ay_{k−1}, k = 1, 2, …  (6)
Alvise Sommariva  Calcolo di autovalori e autovettori  14/ 60

Teorema

Sia A è una matrice quadrata diagonalizzabile avente autovalori λ_{k} tali che

|λ_{1}| > |λ_{2}| ≥ … ≥ |λ_{n}|.

Siano u_{k} ≠ 0 autovettori relativi all’autovalore λ_{k}, cioè

Au_{k} = λ_{k}u_{k}.

Sia

y_{0} = ∑_{k}α_{k}u_{k}, α_{1} ≠ 0.

La successione {y_{s}} definita da y_{s+1} = Ay_{s} converge ad un vettore parallelo a x_{1} e che il coefficiente di Rayleigh (relativo al prodotto scalare euclideo)

ρ(y_{s}, A) := \frac{(y_{s}, Ay_{s})}{(y_{s}, y_{s})} → λ_{1}.  (7)

Dimostrazione.

Per la dimostrazione si confronti [7, p.171]. Essendo la matrice A diagonalizzabile, esistono n autovettori u_{k} (relativi rispettivamente agli autovalori λ_{k}) che sono linearmente indipendenti e quindi formano una base di ℝ^{n}. Sia

y_{0} = ∑_{k}α_{k}u_{k}, α_{1} ≠ 0.

Essendo Au_{k} = λ_{k}u_{k} abbiamo

y_{1} = Ay_{0} = A(∑_{k}α_{k}u_{k}) = ∑_{k}α_{k}Au_{k} = ∑_{k}α_{k}λ_{k}u_{k}

y_{2} = Ay_{1} = A(∑_{k}α_{k}λ_{k}u_{k}) = ∑_{k}α_{k}λ_{k}Au_{k} = ∑_{k}α_{k}λ_{k}^{2}u_{k}

e più in generale y_{s+1} = Ay_{s} = A(∑_{k}α_{k}λ_{k}^{s}u_{k}) = ∑_{k}α_{k}λ_{k}^{s}Au_{k} = ∑_{k}α_{k}λ_{k}^{s+1}u_{k}.

Metodo delle potenze

Osserviamo ora che

\frac{y_{s+1}}{λ_{1}^{s+1}} = ∑_{k}α_{k} \frac{λ_{k}^{s+1}}{λ_{1}^{s+1}} u_{k}  (8)

per cui essendo per k > 1

|\frac{λ_{k}^{s+1}}{λ_{1}^{s+1}} | < 1,

si ha

lim_{s→+∞}(\frac{λ_{k}}{λ_{1}} )^{s} = 0

e quindi la direzione di \frac{y_{s}}{λ_{1}^{s}} , che è la stessa di y_{s}, tende a quella dell’autovettore u_{1}.

Alvise Sommariva  Calcolo di autovalori e autovettori  17/ 60

Si osservi che il coefficiente di Rayleigh ρ(·, A) = (x, Ax)/(x, x) è continuo in ogni x ≠ 0, x ∈ ℝ^{n} in quanto tanto il numeratore quanto il denominatore sono funzioni (multivariate) polinomiali (quadratiche) delle componenti x_{k} di x = (x_{k})_{k} ∈ ℝ^{n}, che sono appunto continue.

Per continuità , se y_{s}/λ^{s} → α_{1}u_{1} allora, essendo λ_{1} ≠ 0, da

lim_{s}ρ(y_{s}, A) := lim_{s} \frac{(y_{s}, Ay_{s})}{(y_{s}, y_{s})} = lim_{s} \frac{(y_{s}/λ_{1}^{s}, A(y_{s}/λ_{1}^{s}))}{(y_{s}/λ_{1}^{s}, y_{s}/λ_{1}^{s})}
= \frac{(α_{1}u_{1}, A(α_{1}u_{1}))}{(α_{1}u_{1}, α_{1}u_{1})} = \frac{(u_{1}, Au_{1})}{(u_{1}u_{1})} = λ_{1}, (9)

ricaviamo che il coefficiente di Rayleigh ρ(y_{s}, A) converge a λ_{1}.

Alvise Sommariva  Calcolo di autovalori e autovettori  18/ 60

Metodo delle potenze

Nota.

Il metodo converge anche nel caso in cui

λ_{1} = … = λ_{r}

per r > 1, tuttavia non è da applicarsi quando l’autovalore di modulo massimo non è unico.

Nota.

In virtú di possibili underflow e underflow si preferisce normalizzare il vettore y_{k} definito in (5). Così l’algoritmo diventa

u_{k} = At_{k−1}  (10)
t_{k} = \frac{u_{k}}{β_{k}} , β_{k} = ∥u_{k}∥_{2}  (11)
l_{k} = ρ(t_{k}, A)  (12)

dove ρ(t_{k}, A) è il coefficiente di Rayleigh definito in (7).

Alvise Sommariva  Calcolo di autovalori e autovettori  19/ 60

Metodo delle potenze inverse

Una variante particolarmente interessante del metodo delle potenze è stata scoperta da Wielandt nel 1944 [9] ed e’ particolarmente utile nel caso in cui A sia una matrice quadrata con n autovettori linearmente indipendenti,

|λ_{1}| ≥ |λ_{2}| ≥ … > |λ_{n}| > 0.  (13)

e si desideri calcolare il piú piccolo autovalore in modulo, cioè λ_{n}, applicando il metodo delle potenze ad A^{−1}.

Metodo delle potenze inverse

Si ottiene cosí la successione {t_{k}} definita da

Au_{k} = t_{k−1}  (14)
t_{k} = \frac{u_{k}}{β_{k}} , β_{k} = ∥u_{k}∥_{2}  (15)

e convergente ad un vettore parallelo a x_{n}. La successione di coefficienti di Rayleigh è tale che

ρ(t_{k}, A^{−1}) := \frac{(t_{k}, A^{−1}t_{k})}{(t_{k}, t_{k})} = \frac{(t_{k}, u_{k+1})}{(t_{k}, t_{k})} → 1/λ_{n}.  (16)

da cui è immediato calcolare λ_{n}.

Alvise Sommariva  Calcolo di autovalori e autovettori  21/ 60

Vediamo in dettaglio questo punto. Se {ξ_{i}} sono gli autovalori di A^{−1} con

|ξ_{1}| > |ξ_{2}| ≥ |ξ_{3}| ≥ … ≥ |ξ_{n}|

allora il metodo delle potenze inverse calcola un’approssimazione di ξ_{1} e di un suo autoversore x.

Si osserva subito che se A^{−1}x = ξ_{i}x (con ξ_{i} ≠ 0) allora moltiplicando ambo membri per ξ_{i}^{−1}A si ottiene leggendo da destra a sinistra Ax = ξ_{i}^{−1}x cioè ξ_{i}^{−1} è un autovalore di A e x è non solo autovettore di A^{−1} relativo all’autovalore ξ_{i}, ma pure autovettore di A relativo all’autovalore ξ_{i}^{−1}.

Conseguentemente se ξ_{1} è l’autovalore di massimo modulo di A^{−1} e λ_{n} è l’autovalore di minimo modulo di A si ha λ_{n} = ξ_{i}^{−1} e che

A^{−1}x = ξ_{1}x ⇒ Ax = ξ_{1}^{−1}x = λ_{n}x

Notiamo che il metodo delle potenze inverse, calcola ξ_{1} = λ_{n}^{−1} e il relativo autovettore x.

Per ottenere λ_{n} viene naturale calcolare ξ_{1}^{−1}, ma usualmente essendo x autovettore di A relativo a λ_{n} si preferisce calcolare λ_{n} via il coefficente di Rayleigh

ρ(x, A) := \frac{(x, Ax)}{(x, x)}.

In generale, fissato μ ∈ ℂ è possibile calcolare, se esiste unico, l’autovalore λ più vicino a μ considerando il seguente pseudocodice [6, p.181]:

(A − μI ) z_{k} = q_{k−1}  (17)
q_{k} = z_{k}/∥z_{k}∥_{2}  (18)
σ_{k} = (q_{k})^{H}Aq_{k}  (19)

Ricordiamo che se λ è autovalore di A allora

Ax = λx ⇒ (A − μI )x = λx − μx = (λ − μ)x

e quindi λ − μ è autovalore di A − μI .

Il metodo delle potenze inverse applicato a A − μI calcola il minimo autovalore σ = λ − μ in modulo di A − μI cioè il σ che rende minimo il valore di |σ| = |λ_{i} − μ|, dove λ_{i} sono gli autovalori di A.

Quindi essendo λ_{i} = σ_{i} − μ si ottiene pure il λ_{i} più vicino a μ. Per versioni piú sofisticate di questa tecnica detta di shift (o in norma infinito invece che in norma 2) si confronti [2, p.379]. Problema. Si può applicare il metodo delle potenze inverse con shift μ nel caso μ sia proprio un autovalore di A?

Metodo QR

Il metodo QR, considerato tra i 10 algoritmi più rilevanti del ventesimo secolo, cerca di calcolare tutti gli autovalori di una matrice A.

Sia A una matrice quadrata di ordine n. Utilizzando il metodo di Gram-Schmidt è possibile fattorizzare la matrice A come prodotto di due matrici Q ed R con Q unitaria (cioè Q^{T} ∗ Q = Q ∗ Q^{T} = I ) ed R triangolare superiore.

Citiamo alcune cose:

▶ La matrice A ha quale sola particolarità di essere quadrata. Nel caso generale però la sua fattorizzazione QR in generale non è unica bensì determinata a meno di una matrice di fase (cf. [2, p.149]) .

▶ Come osservato in [1] p. 614, tale fattorizzazione non è unica (per una opportuna Q unitaria, i segni delle componenti sulla diagonale della matrice R possono essere scelti arbitrariamente). Nel caso sia non singolare, allora tale fattorizzazione è unica qualora si chieda che i coefficienti diagonali di R siano positivi.

▶ La routine Matlab qr effettua tale fattorizzazione. Si consiglia di consultare l’help di Matlab, per consultare le particolarità di tale routine.

Alvise Sommariva  Calcolo di autovalori e autovettori  27/ 60

▶ Se la matrice H è simile a K (cioè esiste una matrice non singolare S tale che H = S^{−1}KS) allora H e K hanno gli stessi autovalori.

Si può vedere facilmente che la relazione di similitudine è transitiva, cioè se H_{1} è simile ad H_{2} e H_{2} è simile ad H_{3} allora H_{1} è simile ad H_{3}.

▶ Dal teorema di Bauer-Fike, segue che se A = V ΛV^{−1} con Λ diagonale, e μ è autovalore di A + Δ, allora esiste λ autovalore di A tale che |λ − μ| ≤ K_{∞}(A)∥Δ∥_{∞}.

Il metodo QR venne pubblicato indipendemente nel 1961 da Francis e da Kublanovskaya e successivamente implementato in EISPACK. Ci limiteremo a considerare versioni di base del metodo.

Sia

A_{0} = A = Q_{0}R_{0}

e si ponga

A_{1} := R_{0}Q_{0}.

Poichè

Q_{0}A_{1}Q_{0}^{T} = Q_{0}A_{1}Q_{0}^{T} = Q_{0}R_{0}Q_{0}Q_{0}^{T} = A_{0}

la matrice A_{1} è simile ad A_{0} (si ponga S = Q_{0}^{−1} = Q_{0}^{T} ) e quindi ha gli stessi autovalori. Sia quindi in generale

A_{k} = Q_{k}R_{k}

A_{k+1} = R_{k}Q_{k}.

Per le stesse motivazioni A_{k+1} è simile ad A_{k}, e per transitività ad A_{0}. Quindi A_{k+1} ha gli stessi autovalori di A_{0}.

Alvise Sommariva  Calcolo di autovalori e autovettori  29/ 60

Metodo QR

Per la convergenza del metodo esistono vari risultati (cf. [3, p.393], [4, p.352], [7, p.180]). Da [6, p.169]

Teorema

Se A ∈ ℝ^{n×n} ha autovalori tutti distinti in modulo, con

|λ_{1}| > … > |λ_{n}|  (20)

allora l’alg. QR converge a A_{∞} = (a_{i,j}^{∞}) triangolare sup., cioè

lim_{k}A_{k} = ( \begin{matrix}a_{1,1}^{∞} & a_{1,2}^{∞} & … & … & a_{1,n}^{∞} \\ 0 & a_{2,2}^{∞} & a_{2,3}^{∞} & … & a_{2,n}^{∞} \\ 0 & 0 & a_{3,3}^{∞} & … & a_{3,n}^{∞} \\ 0 & 0 & … & … & … \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & 0 & 0 & a_{n−1,n−1}^{∞} & a_{n−1,n}^{∞} \\ 0 & 0 & 0 & 0 & a_{n,n}^{∞} \end{matrix} )  (21)

con λ_{k} = a_{k,k}^{∞}.

Alvise Sommariva  Calcolo di autovalori e autovettori  30/ 60

Inoltre

▶ Se la condizione (20) non è verificata si può dimostrare che la successione {A_{k}} tende a una forma triangolare a blocchi.

▶ se A_{k} = (a_{i,j}^{(k)}), e λ_{i−1} ≠ 0

|a_{i,i−1}^{(k)}| = O (\frac{|λ_{i}|}{|λ_{i−1}|})^{k} , i = 2, …, n, k → ∞.  (22)

▶ se la matrice è simmetrica, allora

A_{∞} = diag(λ_{1}, …, λ_{1}).

▶ se A è una matrice Hessenberg superiore allora l’algoritmo QR converge ad A_{∞} triangolare a blocchi, simile ad A e con gli autovalori di ogni blocco diagonale tutti uguali in modulo.

Alvise Sommariva  Calcolo di autovalori e autovettori  31/ 60

Implementazione del metodo QR

Nelle implementazioni si calcola con un metodo scoperto da Householder (ma esiste un metodo alternativo dovuto a Givens) una matrice di Hessenberg T

T = ( \begin{matrix}a_{1,1} & a_{1,2} & a_{1,3} & … & a_{1,n} \\ a_{2,1} & a_{2,2} & a_{2,3} & … & a_{2,n} \\ 0 & a_{3,2} & a_{3,3} & … & a_{3,n} \\ 0 & 0 & a_{4,3} & … & a_{4,n} \\ … & … & … & … & … \\ 0 & 0 & 0 & a_{n,n−1} & a_{n,n} \end{matrix} )

simile ad A ed in seguito si applica il metodo QR relativamente alla matrice T . Se A è simmetrica la matrice T risulta tridiagonale simmetrica.

Si può mostrare che se A è una matrice di Hessenberg superiore, allora A = QR con Q di Hessenberg superiore.

Inoltre, se A è tridiagonale allora A = QR con Q di Hessenberg e R triangolare superiore con r_{i,j} = 0 qualora j − i ≥ 2.

In entrambi i casi le iterazioni mantengono la struttura, cioè se A_{0} = T è di Hessenberg, allora A_{k} è di Hessenberg, se A_{0} = T è tridiagonale allora A_{k} è tridiagonale.

Il numero di moltiplicazioni necessarie

▶ all’algoritmo di Givens per calcolare tale matrice T a partire

da A è approssimativamente 10n^{3}/3; ▶ all’algoritmo di Householder per calcolare tale matrice T a

partire da A è approssimativamente 5n^{3}/3. Il metodo QR applicato ad una matrice A in forma di Hessenberg superiore ha ad ogni passo un costo di 2n^{2} operazioni moltiplicative.

Per versioni piú sofisticate come il metodo QR con shift, si veda [3], p. 394.

Il metodo delle potenze in Matlab

Partiamo con una versione semplice power basic del metodo delle potenze

 f u n c t i o n [ lambda1 , x1 , niter , err]= power_basic ( A , z0 , toll , nmax )

% INPUT :
% A  : MATRICE DI CUI VOGLIAMO CALCOLARE L ’AUTOVALORE DI MASSIMO MODULO.
% z0 : VETTORE INIZIALE (NON NULLO) .
% t o l l : TOLLERANZA .
% nmax : NUMERO MASSIMO DI ITERAZIONI .

%

% OUTPUT:
% lambda1 : VETTORE DELLE APPROSSIMAZIONI DELL ’AUTOVALORE DI MASSIMO MODULO.
% x1  : AUTOVETTORE RELATIVO ALL ’AUTOVALORE DI MASSIMO MODULO.
% n i t e r  : NUMERO DI ITERAZIONI .
% e r r  : VETTORE DEI RESIDUI PESATI RELATIVI A ” lambda1 ” .

%

% TRATTO DA QUARTERONI−SALERI , ”MATEMATICA NUMERICA” , p . 1 8 4 .

%

Il metodo delle potenze in Matlab

q=z0/norm ( z0 ) ; q2=q ; err = [ ] ; lambda1 = [ ] ;

res=toll +1; niter =0; z=A∗q ;
w h i l e ( res >= toll & niter <= nmax )

q=z/norm ( z ) ; z=A∗q ; lam=q ’ ∗ z ; x1=q ;

z2=q2 ’ ∗ A ; q2=z2/norm ( z2 ) ; q2=q2 ’ ; y1=q2 ; costheta=abs ( y1 ’ ∗ x1 ) ;
niter=niter +1; res=norm ( z−lam∗q ) / costheta ;
err=[err ; res ] ; lambda1 =[ lambda1 ; lam ] ;

end

Qualche nota

▶ il vettore iniziale z0 e’ normalizzato ed in err, lambda1 vengono memorizzati rispettivamente i valori dell’errore compiuto e dell’autovalore di massimo modulo λmax;

▶ l’assegnazione res=toll+1; forza l’algoritmo ad entrare nel ciclo while, mentre z=A*q; è una quantità da utilizzarsi per il calcolo dell’autovalore λmax; ▶ nel ciclo while, q è un’approssimazione di un autoversore

relativo a λmax, mentre lam di λmax;

▶ il ciclo si interrompe se un numero massimo di iterazioni niter è raggiunto oppure

\frac{||Aq^{k} − λ^{k}||_{2}}{| cos(θ_{λ_{k}})|} < tol

dove θ_{λ_{k}} è l’angolo formato tra (un’approssimazione del)l’autovalore destro x1 e sinistro y1 associati a lam (cf. [6, p.180])

Il metodo delle potenze in Matlab: esempio I

Testiamo il codice per il calcolo dell’autoval. di massimo modulo di

A = ( \begin{matrix}−15.5 & 7.5 & 1.5 \\ −51 & 25 & 3 \\ −25.5 & 7.5 & 11.5 \end{matrix} )

= ( \begin{matrix}1 & 2 & 3 \\ 2 & 5 & 6 \\ 7 & 9 & 3 \end{matrix} ) · ( \begin{matrix}10 & 0 & 0 \\ 0 & 10 & 0 \\ 0 & 0 & 1 \end{matrix} ) · ( \begin{matrix}1 & 2 & 3 \\ 2 & 5 & 6 \\ 7 & 9 & 3 \end{matrix} )^{−1}(23)

La matrice A è diagonalizzabile e ha autovalori 10, 10, 1. Si può vedere che una base di autovettori relativa agli autovalori 10, 10, 1 è composta da (1, 2, 7), (2, 5, 9), (3, 6, 3). Quale vettore iniziale del metodo delle potenze consideriamo

z_{0} = (1, 1, 1) = (7/6) · (1, 2, 7) − 1 · (2, 5, 9) + (11/18) · (3, 6, 3)

e quindi il metodo delle potenze applicato ad A, e avente quale punto iniziale z_{0} può essere utilizzato per il calcolo dell’autovalore di massimo modulo di A, poichè α_{1} = 7/6 ≠ 0.

Alvise Sommariva  Calcolo di autovalori e autovettori  39/ 60

Il metodo delle potenze in Matlab: esempio I

Dalla shell di Matlab/Octave:

>> S=[1 2 3 ; 2 5 6 ; 7 9 3 ] ;
>> D=d i a g ( [ 1 0 10 1 ] ) ;
>> A=S∗D∗ i n v ( S )

 A =

−15.5000  7.5000  1.5000
−51.0000  25.0000  3.0000
−25.5000  7.5000  11.5000
>> z0 =[1 1 1 ] ’ ;
>> toll=10ˆ (−8) ;
>> nmax =10;
>> format short e ;

>> [ lambda1 , x1 , niter , err]= power_basic ( A , z0 , toll , nmax )
lambda1 =

1.1587 e+001
1.0138 e+001
1.0014 e+001
1.0001 e+001

…

1.0000 e+001
1.0000 e+001
1.0000 e+001

x1 =

−2.8583e−001
−9.1466e−001
−2.8583e−001

niter =

10

err =

2.2466 e+000
2.1028 e−001
2.0934 e−002
2.0925 e−003

?

2.0924 e−007
2.0924 e−008
2.0924 e−009

>>

La convergenza è abbastanza veloce come si vede dalla quantità err, che consiste in un particolare residuo pesato. Una questione sorge spontanea. Cosa sarebbe successo se avessimo utilizzato l’algoritmo senza normalizzazione come ad esempio

power method definito da

f u n c t i o n [ lambda , v]= power_method ( A , x0 , maxit )
v=x0 ;
f o r index =1: maxit
     v_old=v ;
     v=A∗ v_old ;
     lambda=(v_old ’ ∗ v ) /( v_old ’ ∗ v_old ) ;

end

Proviamo il test, facendo iterare il metodo prima 5, poi 100 volte e alla fine 1000 volte (si noti il settaggio della variabile maxit relativa al numero di iterazioni da compiere):

>> x0 =[1 1 1 ] ’

x0 =

1 1 1

>> A=[−15.5 7 . 5 1 . 5 ; −51 25 3 ; −25.5 7 . 5 1 1 . 5 ]

A =

−15.5000  7.5000  1.5000
−51.0000  25.0000  3.0000
−25.5000  7.5000  11.5000

>> [ lambda , v]= power_method ( A , x0 , 5 )
lambda =

10.0014

v =

1 . 0 e+005 ∗

−0.8333 −2.6666 −0.8333

>> [ lambda , v]= power_method ( A , x0 , 1 0 0 )
lambda =

10.0000

v =

1 . 0 e+100 ∗

−0.8333 −2.6667 −0.8333

>> [ lambda , v]= power_method ( A , x0 , 1 0 0 0 )
lambda =

NaN

v =

NaN NaN NaN

>>

La ragione è semplice. Per k relativamente piccolo si ha A · t_{k} ≈ 10 · t_{k} e quindi per s ≥ k

t_{s} ≈ A^{s−k} · t_{k} ≈ 10 · A^{s−k−1} · t_{k} ≈ ⋯ ≈ 10^{s−k} · t_{k}

da cui

∥t_{s}∥_{2} ≈ 10^{s−k} · ∥t_{k}∥_{2}

spiegando quindi perchè si possano avere problemi di overflow applicando l’algoritmo di base.

Il metodo delle potenze in Matlab: esempio II

Proviamo un test diverso, questa volta con la matrice (diagonalizzabile)

A = ( \begin{matrix}1 & 2 \\ 0 & −1 \end{matrix} ) ,

avente autovalori λ_{1} = 1 e λ_{2} = −1 e autovettori linearmente indipendenti (1, 0), (−1, 1). Quale vettore iniziale poniamo

x_{0} = (1, 3) = 4 · (1, 0) + 3 · (−1, 1)

e quindi il metodo delle potenze applicato ad A, partendo da x_{0} può essere sicuramente applicato. D’altra parte dubitiamo converga in quanto |λ_{1}| = |λ_{2}| = 1 pur essendo λ_{1} ≠ λ_{2}.

Il metodo delle potenze in Matlab: esempio II

Dalla shell di Matlab/Octave:

>> A=[1 2 ; 0 −1]

A =

1  2
0  −1

>> [ lambda1 , x1 , niter , err]= power_basic ( A , [ 1 ; 3] ,10ˆ ( −8) ,15)
lambda1 =

−3.4483e−002
−2.0000e−001

…

−3.4483e−002
−2.0000e−001
−3.4483e−002
−2.0000e−001

x1 =

3.1623 e−001
9.4868 e−001

niter =

16

err =

4.4567 e−001
2.4000 e+000

…

4.4567 e−001
2.4000 e+000
4.4567 e−001
2.4000 e+000

>>

Dal residuo pesato è chiaro che il metodo non converge, e come anticipato il motivo è la presenza di autovalori distinti aventi modulo massimo.

Alvise Sommariva  Calcolo di autovalori e autovettori  47/ 60

Il metodo delle potenze in Matlab: esempio III

Vediamo il caso della matrice diagonalizzabile (autovalori distinti!)

A = ( \begin{matrix}1 & 2 \\ 0 & 10 \end{matrix} ) ,

in cui il metodo funziona rapidamente (λmax = 10).

>> A=[1 2 ; 0 1 0 ] ; [ lambda1 , x1 , niter , err]= power_basic ( A , [ 1 ; 3] ,10ˆ ( −8) ,15)

lambda1 =

9.9779 e+000
9.9979 e+000
9.9998 e+000
1.0000 e+001
1.0000 e+001
1.0000 e+001
1.0000 e+001
1.0000 e+001

x1 =

2.1693 e−001
9.7619 e−001

niter =

8

err =

9.6726 e−002
9.7529 e−003
9.7610 e−004
9.7618 e−005
9.7619 e−006
9.7619 e−007
9.7619 e−008
9.7619 e−009

Si osservi che il metodo termina in quanto l’errore pesato err è minore della tolleranza toll = 10^{−8}.

Alvise Sommariva  Calcolo di autovalori e autovettori  48/ 60

Il metodo delle potenze inverse in Matlab

Una versione di base invpower del metodo delle potenze inverse [6, p.184] è

 f u n c t i o n [ lambda , x , niter , err]= invpower ( A , z0 , mu , toll , nmax )

% DATO UN VALORE mu, SI CALCOLA L ’AUTOVALORE ”lambda mu” PIU ’ VICINO A mu .
% INPUT :
% A  : MATRICE DI CUI VOGLIAMO CALCOLARE L ’AUTOVALORE ”lambda mu ” .
% z0 : VETTORE INIZIALE (NON NULLO) .
% mu : VALORE DI CUI VOGLIAMO CALCOLARE L ’AUTOVALORE PIU ’ VICINO .
% t o l l : TOLLERANZA .
% nmax : NUMERO MASSIMO DI ITERAZIONI .

%

% OUTPUT:
% lambda : VETTORE DELLE APPROSSIMAZIONI DELL ’AUTOVALORE DI MINIMO MODULO.
% x  : AUTOVETTORE RELATIVO ALL ’AUTOVALORE DI MINIMO MODULO.
% n i t e r : NUMERO DI ITERAZIONI .
% e r r  : VETTORE DEI RESIDUI PESATI RELATIVI A ” lambda ” .

%

% TRATTO DA QUARTERONI−SALERI , ”MATEMATICA NUMERICA” , p . 1 8 4 .

%

Il metodo delle potenze inverse in Matlab

n=max( s i z e ( A ) ) ; M=A−mu∗ eye ( n ) ; [ L , U , P]= l u ( M ) ;
q=z0/norm ( z0 ) ; q2=q ’ ; err = [ ] ; lambda = [ ] ;

res=toll +1; niter =0;
w h i l e ( res >= toll & niter <= nmax )
niter=niter +1; b=P∗q ; y=L\b ; z=U\y ;

q=z/norm ( z ) ; z=A∗q ; lam=q ’ ∗ z ;
b=q2 ’ ; y=U ’\ b ; w=L ’\ y ;
q2=(P ’ ∗ w ) ’ ; q2=q2/norm ( q2 ) ; costheta=abs ( q2∗q ) ;

i f ( costheta > 5e−2)
     res=norm ( z−lam∗q ) / costheta ; err=[err ; res ] ; lambda =[ lambda ; lam ] ;

e l s e

d i s p ( ’\n \t [ ATTENZIONE ] : AUTOVALORE MULTIPLO ’ ) ; break ;

end

x=q ;

end

Il metodo delle potenze inverse in Matlab

Forniamo ora alcune spiegazioni del codice in invpower.

▶ Per risolvere il sistema lineare in 17, si effettua una fattorizzazione PM = LU della matrice M = A − μI ;

▶ All’interno del ciclo while, nella prima riga si calcola z_{k}, mentre nella successiva un suo versore q_{k}, e σ_{k} è immagazzinato in lam;

▶ Similmente al metodo diretto si effettua il prodotto scalare di un’autovalore sinistro con uno destro.

Il metodo delle potenze inverse in Matlab: esempio I

Applichiamo il metodo delle potenze inverse per il calcolo dell’autovalore più piccolo in modulo della matrice

A = ( \begin{matrix}−15.5 & 7.5 & 1.5 \\ −51 & 25 & 3 \\ −25.5 & 7.5 & 11.5 \end{matrix} )

= ( \begin{matrix}1 & 2 & 3 \\ 2 & 5 & 6 \\ 7 & 9 & 3 \end{matrix} ) · ( \begin{matrix}10 & 0 & 0 \\ 0 & 10 & 0 \\ 0 & 0 & 1 \end{matrix} ) · ( \begin{matrix}1 & 2 & 3 \\ 2 & 5 & 6 \\ 7 & 9 & 3 \end{matrix} )^{−1}(24)

Come visto la matrice A è quindi diagonalizzabile, ha autovalori 10, 10, 1 e relativi autovettori è (1, 2, 7), (2, 5, 9), (3, 6, 3) formanti una base di ℝ^{3}. Quale vettore iniziale consideriamo

z_{0} = (1, 1, 1) = (7/6) · (1, 2, 7) − 1 · (2, 5, 9) + (11/18) · (3, 6, 3)

e quindi il metodo delle potenze inv. applicato ad A, e avente quale punto iniziale z_{0} può essere utilizzato per il calcolo dell’autovalore di minimo modulo di A.

Alvise Sommariva  Calcolo di autovalori e autovettori  52/ 60

>> z0 = [ 1 ; 1 ; 1 ] ; mu =0; toll=10ˆ (−8) ; nmax =10;
>> A=[−15.5 7 . 5 1 . 5 ; −51 25 3 ; −25.5 7 . 5 1 1 . 5 ] ;

>> [ lambda , x , niter , err]= invpower ( A , z0 , mu , toll , nmax )
lambda =

0.39016115351993 0.94237563941268 0.99426922936854 0.99942723776656 0.99994272692315 0.99999427272378 0.99999942727270 0.99999994272728 0.99999999427273

x =

0.40824829053809 0.81649658085350 0.40824829053809

niter =

9

err =

0.81535216507377 0.08358101289062 0.00838126258396 0.00083836078891 0.00008383842712 0.00000838386620 0.00000083838685 0.00000008383868 0.00000000838387

>>

La convergenza è lineare (come si intuisce dalle approssimazioni contenute nel vettore lambda).

Alvise Sommariva  Calcolo di autovalori e autovettori  53/ 60

Il metodo delle potenze inverse in Matlab: esempio I

Figura: Grafico che illustra la convergenza lineare del metodo delle potenze inverse nell’esempio 1.

Alvise Sommariva  Calcolo di autovalori e autovettori  54/ 60

Per vederlo, dalla shell di Matlab/Octave calcoliamo l’errore assoluto/relativo relativo all’autovalore 1:

>> s=1−lambda
s =

0.60983884648007 0.05762436058732 0.00573077063146 0.00057276223344 0.00005727307685 0.00000572727622 0.00000057272730 0.00000005727272 0.00000000572727

>> s e m i l o g y ( 1 : l e n g t h ( s ) , s )

generando il grafico in scala semi-logaritmica in figura che evidentemente sottolinea la convergenza lineare.

Il metodo QR in Matlab

Una versione di base del metodo QR è la seguente. Si salvi il files houshess.m che implementa la trasformazione per similitudine di A in una matrice di Hessenberg

f u n c t i o n [ H , Q]= houshess ( A )
% REDUCTION OF A MATRIX TO A SIMILAR HESSENBERG ONE.
% SEE QUARTERONI, SACCO, SALERI P . 1 9 2 .
n=max( s i z e ( A ) ) ; Q=eye ( n ) ; H=A ;
f o r k =1:( n−2)
[ v , beta ]= vhouse ( H ( k+1:n , k ) ) ; I=eye ( k ) ; N=z e r o s ( k , n−k ) ;
m=l e n g t h ( v ) ; R=eye ( m )−beta ∗v∗v ’ ; H ( k+1:n , k : n )=R∗H ( k+1:n , k : n ) ;
H ( 1 : n , k+1:n )=H ( 1 : n , k+1:n ) ∗R ; P=[I , N ; N ’ , R ] ; Q=Q∗P ;

end

Il metodo QR in Matlab

ove vhouse.m è definito da

f u n c t i o n [ v , beta ]= vhouse ( x )
% BUILDING HOUSEHOLDER VECTOR.
% SEE QUARTERONI, SACCO, SALERI P . 1 9 7 .
n=l e n g t h ( x ) ; x=x/norm ( x ) ; s=x ( 2 : n ) ’∗ x ( 2 : n ) ; v =[1; x ( 2 : n ) ] ;
i f ( s==0)
beta =0;

e l s e

mu=s q r t ( x ( 1 ) ˆ 2+s ) ;
i f ( x ( 1 ) <= 0)
v ( 1 )=x ( 1 )−mu ;

e l s e

v ( 1 )=−s /( x ( 1 )+mu ) ;

end

beta=2∗v ( 1 ) ˆ 2/( s+v ( 1 ) ˆ 2) ;
v=v/v ( 1 ) ;

end

Il metodo QR in Matlab

Quindi QRbasicmethod.m che calcola la matrice triangolare T relativa a QR

f u n c t i o n [ T , h i s t ]= QRbasicmethod ( T_input , maxit )

% QR METHOD FOR A SYMMETRIC TRIDIAGONAL MATRIX ” T i n p u t ” .

T=T_input ;

h i s t=s o r t ( d i a g ( T ) ) ;

f o r index =1: maxit

[ Q , R]= qr ( T ) ;
T=R∗Q ;  % NEW SIMILAR MATRIX .
h i s t =[ h i s t s o r t ( d i a g ( T ) ) ] ; % IT STORES THE DIAGONAL ELEMENTS
% OF THE ” i n d e x ” ITERATION .

end

Esercizio

▶ Data la matrice di Hilbert di ordine 5, ottenibile in Matlab col comando hilb(5) si calcolino col metodo delle potenze i suoi minimi e massimi autovalori in modulo.

▶ Da questi si determini il condizionamento della matrice in norma 2 e lo si confronti con cond(hilb(5),2). Eseguire lo stesso esercizio utilizzando il metodo QR.

Bibliografia

K. Atkinson, Introduction to Numerical Analysis, Wiley, 1989. D. Bini, M. Capovani e O. Menchi, Metodi numerici per l’algebra lineare, Zanichelli, 1988. V. Comincioli, Analisi Numerica, metodi modelli applicazioni, Mc Graw-Hill, 1990. G.H. Golub e C.F. Van Loan, Matrix Computation, 3rd Edition, The John Hopkins University Press 1996. Netlib,http://www.netlib.org/templates/matlab/ A. Quarteroni, R. Sacco, F. Saleri, Matematica numerica, 2001. A. Quarteroni e F. Saleri, Introduzione al calcolo scientifico, Springer Verlag, 2006. Wikipedia (Matrice irriducibile) http://it.wikipedia.org/wiki/Matrice irriducibile Wikipedia (Inverse iteration) http://en.wikipedia.org/wiki/Inverse iteration. Wikipedia (Metodo delle potenze) http://it.wikipedia.org/wiki/Metodo delle potenze. Wikipedia (PageRank) http://it.wikipedia.org/wiki/Page rank. Wikipedia (Rayleigh quotient) http://en.wikipedia.org/wiki/Rayleigh quotient. Wikipedia (QR algorithm) http://en.wikipedia.org/wiki/QR algorithm. Wikipedia (QR decomposition) http://en.wikipedia.org/wiki/QR decomposition. Wikipedia (Teoremi di Gershgorin) http://it.wikipedia.org/wiki/Teoremi di Gershgorin.

Alvise Sommariva  Calcolo di autovalori e autovettori  60/ 60