Autovalori

20 febbraio 2009

1 Calcolo degli autovalori di una matrice

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

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) mentre in altri casi si vogliono determinare tutti gli autovalori.

Per semplicità , dopo i teoremi di localizzazione di Gerschgorin, 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 [3].

Una interessante applicazione è l’algoritmo di PageRank [16], 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 (cf. [1], [4], [5]).

2 Teoremi di Gerschgorin

In questo paragrafo mostriamo tre teoremi di localizzazione di autovalori dovuti a Gerschgorin (cf. [3, p.76], [13].

Teorema 2.1 (Primo teorema di Gerschgorin). Gli autovalori di una matrice di ordine sono tutti contenuti nell’unione dei cerchi di Gerschgorin

� : _{}_{1}^{}_{}

Vediamo quale esempio la matrice

{ \atop }\begin{matrix}15 & 2 & 2 \\ 1 & 10 & 3 \\ 2 & 1 & 0 \end{matrix} { \atop }(2)

\underline{Autovalori}

Figura 1: Cerchi di Gerschgorin della matrice A definita in (2)

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

1  � : 15  2 2 4
2  � : 10 1  3 4
3  � : 0  2 1 3

Teorema 2.2 (Secondo teorema di Gerschgorin). Se l’unione _{1} di cerchi di Gerschgorin è disgiunta dall’unione _{2} dei rimanenti , allora autovalori appartengono a _{1} e appartengono a _{2}.

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

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

Questo teorema sarà utile in seguito per dimostrare che una matrice, detta di Poisson, è non singolare. Ricordiamo che una matrice si dice irriducibile se non è riducibile e che una matrice di ordine 2 è riducibile se esiste una matrice di permutazione e un intero , 0 , tale che

\begin{matrix}_{1 1} & _{1 2} \\ 0 & _{2 2} \end{matrix}

in cui _{1 1} �^{}, _{12 2} �^{( ) ( )} .

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

>> eig(A)
ans =
    0.5121
   14.1026
   10.3854
>>

a conferma di quanto stabilito dai primi due teoremi di Gerschgorin. Nota. Cosa possiamo dire relativamente agli autovalori di se applichiamo i teoremi di Gerschgorin ad ^{}invece che ad ?

3 Il 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 una matrice quadrata di ordine con autovettori _{1}, , _{}linearmente indipendenti e autovalori _{1}, , _{}tali che

1  _{2}  _{}(3)

A tal proposito ricordiamo (cf. [6], p. 951) che

1. una matrice d́iagonalizzabile se e solo se possiede autovettori linearmente indipendenti;

2. Se tutti gli autovalori di sono distinti la matrice è diagonalizzabile;

3. Una matrice simmetrica (hermitiana) è diagonalizzabile.

Il metodo delle potenze è definito come segue. Sia _{0} ^{}definito da

0 _{1}^{}_{}_{}_{1} 0

e si generi la successione

0  _{0}  (4)
1  1 2  (5)

Teorema 3.1 Sia è una matrice quadrata diagonalizzabile avente autovalori _{} tali che

1  _{2}  _{}

Siano _{}0 autovettori relativi all’autovalore _{}, cioè Sia

0 _{}_{}_{}_{1} 0

La successione _{}definita da _{1} _{}converge ad un vettore parallelo a _{1} e che il coefficiente di Rayleigh (relativo al prodotto scalare euclideo)

( _{}) : \frac{( _{}_{})}{( _{}_{})}  (6)

converge a _{1}.

Dimostrazione (cf. [11, p.171]). Essendo la matrice diagonalizzabile, esistono autovettori _{}(relativi rispettivamente agli autovalori _{}) che sono linearmente indipendenti e quindi formano una base di �^{}. Sia

0 _{}_{}_{}_{1} 0

Essendo _{}_{}_{}abbiamo

1 _{0} (_{}_{}_{}) _{}_{}_{}_{}_{}_{}_{}

2 _{1} (_{}_{}_{}_{}) _{}_{}_{}_{}_{}_{}_{}^{2} _{}

e più in generale

1 _{}(_{}_{}_{}^{}_{}) _{}_{}_{}^{}_{}_{}_{}_{}^{1} _{}

Osserviamo ora che

\frac{_{1}}{_{1}^{}}_{}_{}\frac{_{}^{1}}{_{1}^{}}_{}(7)

per cui essendo

; ; ; ; \frac{_{}^{1}}{_{1}^{}} ; ; ; ; 1 1

la direzione di \frac{_{1}}{_{1}^{}}tende a quella di _{}. Per continuità il coefficiente di Rayleigh

( _{}) : \frac{( _{}_{})}{( _{}_{})}  (8)

converge a _{1}.

Si osservi che il coefficiente di Rayleigh ( ) è continuo in ogni 0, �^{} in quanto tanto il numeratore quanto il denominatore sono funzioni polinomiali (quadratiche) delle componenti _{}di ( _{})_{}�^{}, che sono appunto continue.

Il metodo converge anche nel caso in cui

1 _{} per 1, tuttavia non è da applicarsi quando l’autovalore di modulo massimo non è unico. In virtú di alcune possibili problemi di underflow e underflow si preferisce normalizzare passo passo il vettore _{}definito in (4). Conseguentemente l’algoritmo del metodo delle potenze risulta

1  (9)
\frac{_{}}{_{}}_{}_{}_{2}  (10)
( _{})  (11)

dove ( _{}) è il coefficiente di Rayleigh definito in (6).

4 Il metodo delle potenze inverse

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

1  _{2}  _{}0  (12)

e si desideri calcolare il piú piccolo autovalore in modulo, cioè _{}, applicando il metodo delle potenze ad ^{1} . Si ottiene cosí la successione _{}definita da

1  (13)
\frac{_{}}{_{}}_{}_{}_{2}  (14)

e convergente ad un vettore parallelo a _{}. La successione di coefficienti di Rayleigh

( _{}^{1} ) : \frac{( _{}^{1} _{})}{( _{}_{})} \frac{( _{}_{1} )}{( _{}_{})} 1/ _{}(15)

da cui è immediato calcolare _{}.

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

1  _{2}  _{3}  _{}

allora il metodo delle potenze inverse calcola un’approssimazione di _{1} e di un suo autoversore . Si osserva subito che se ^{1} _{}(con _{}0) allora moltiplicando ambo membri per _{}^{1} si ottiene leggendo da destra a sinistra _{}^{1} cioè _{}^{1} è un autovalore di e è non solo autovettore di ^{1} relativo all’autovalore , ma pure autovettore di relativo all’autovalore _{}^{1} . Conseguentemente se _{1} è l’autovalore di massimo modulo di ^{1} e _{}è l’autovalore di minimo modulo di si ha _{}_{}^{1} e che

1  _{1}  _{1}^{1}  _{}

Notiamo che il metodo delle potenze inverse, calcola _{1} _{}^{1} e il relativo autovettore . Per ottenere _{}viene naturale calcolare _{1}^{1} , ma usualmente essendo autovettore di relativo a _{}si preferisce per questioni numeriche calcolare _{}via coefficente di Rayleigh

( ) : \frac{( )}{( )}

In generale, fissato � è possibile calcolare l’autovalore più vicino a considerando il seguente pseudocodice [10, p.181]:

(  ) _{}_{1}  (16)
/ _{}_{2}  (17)
( _{})^{}_{}(18)

Ricordiamo che se è autovalore di allora

(  )  ( )
e quindi  è autovalore di  . Il metodo delle potenze inverse applicato a
calcola il minimo autovalore  in modulo di  cioè il che rende

minimo il valore di _{}, dove _{}sono gli autovalori di . Quindi essendo si ottiene pure il _{}più vicino a .

Per versioni piú sofisticate di questa tecnica detta di shift (o in norma infinito invece che in norma 2) si confronti [3, p.379]. Problema. Si può applicare il metodo delle potenze inverse con shift nel caso sia proprio un autovalore di A?

5 Il metodo QR

Sia una matrice quadrata di ordine . Utilizzando il metodo di Householder è possibile fattorizzare la matrice A come prodotto di due matrici ed con unitaria (cioè ^{}^{}) ed triangolare superiore.

Citiamo alcune cose:

1. La matrice ha quale sola particolarità di essere quadrata.

2. Come osservato in [2] p. 614, tale fattorizzazione non è unica (i segni delle componenti sulla diagonale della matrice possono essere scelti arbitrariamente).

3. La routine Matlab qr effettua tale fattorizzazione.

4. Se la matrice è simile a (cioè esiste una matrice non singolare tale che 1 ) allora e hanno gli stessi autovalori. Si può vedere facilmente che la relazione di similitudine è transitiva, cioè se _{1} è simile ad _{2} e _{2} è simile ad _{3} allora _{1} è simile ad _{3}.

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

Sia

0 _{0} _{0}

e si ponga

1 : _{0} _{0}

Poichè

0 _{1} _{0}^{}_{0} _{1} _{0}^{}_{0} _{0} _{0} _{0}^{}_{0}

la matrice _{1} è simile ad _{0} (si ponga _{0}^{1} _{0}^{}) e quindi ha gli stessi autovalori. Sia quindi in generale

1 _{}_{}

Per le stesse motivazioni _{1} è simile ad _{}, e per transitività ad _{0}. Quindi _{1} ha gli stessi autovalori di _{0}. Per la convergenza del metodo citiamo il seguente risultato (cf. [6], p. 393).

Teorema 5.1 Sia ^{}la matrice formata dagli autovettori di sinistra della matrice . Se tutti i minori principali di sono diversi da 0 e se gli autovettori di verificano la condizione

1  _{2}  _{}0  (19)

allora la matrice _{1} _{}^{}_{}tende a una forma triangolare superiore e ( _{1} )_{} converge verso _{}.

Se la condizione (19) non è verificata si può dimostrare che la successione _{}tende a una forma triangolare a blocchi. Per versioni piú sofisticate come il metodo QR con shift, si veda [6], p. 394. Nelle implementazioni si calcola con un metodo scoperto da Householder (ma esiste un metodo alternativo dovuto a Givens) una matrice di Hessenberg

; ; ; ; ; ; ; ; \begin{matrix}_{1 1} & _{1 2} & _{1 3} & & _{1} \\ _{2 1} & _{2 2} & _{2 3} & & _{2} \\ 0 & _{3 2} & _{3 3} & & _{3} \\ 0 & 0 & _{4 3} & & _{4} \\ & & & & \\ 0 & 0 & 0 & _{1} & _{} \end{matrix} ; ; ; ; ; ; ; ;

simile ad ed in seguito si applica il metodo QR relativamente alla matrice . Se è simmetrica la matrice risulta tridiagonale simmetrica. Per inciso, si può dimostrare che la velocità di convergenza dipende dal rapporto

max_{1 1} max \frac{_{1}}{_{}}

6 Esperimenti numerici in Matlab

6.1 Il metodo delle potenze

Partiamo con una versione semplice power_basic del metodo delle potenze

function [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).
% toll: TOLLERANZA.
% nmax: NUMERO MASSIMO DI ITERAZIONI.
%
% OUTPUT:
% lambda1 : VETTORE DELLE APPROSSIMAZIONI DELL'AUTOVALORE DI MASSIMO MODULO.
% x1      : AUTOVETTORE RELATIVO ALL'AUTOVALORE DI MASSIMO MODULO.
% niter   : NUMERO DI ITERAZIONI.
% err     : VETTORE DEI RESIDUI PESATI RELATIVI A "lambda1".
%
% TRATTO DA QUARTERONI-SALERI, "MATEMATICA NUMERICA", p. 184.
%
q=z0/norm(z0); q2=q; err=[]; lambda1=[];
res=toll+1; niter=0; z=A*q;
while (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

1. 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;

2. 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;

3. nel ciclo while, q è un’approssimazione di un autoversore relativo a max, mentre lam di max;

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

\frac{^{}^{}_{2}}{cos( _{_{}})} tol

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

6.2 Esempio 1

Testiamo il codice relativamente al calcolo dell’autovalore di massimo modulo di

{ \atop }\begin{matrix}15 5 & 7 5 & 1 5 \\ 51 & 25 & 3 \\ 25 5 & 7 5 & 11 5 \end{matrix} { \atop }

{ \atop }1 2 3 2 5 6 7 9 3 { \atop }{ \atop }10 0 0 0 10 0 0 0 1 { \atop }{ \atop }1 2 3 2 5 6 7 9 3 { \atop }^{1}  (20)

La matrice è 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

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 , e avente quale punto iniziale _{0} può essere utilizzato per il calcolo dell’autovalore di massimo modulo di , poichè 1 7/6 0. Dalla shell di Matlab/Octave:

>> S=[1 2 3; 2 5 6; 7 9 3];
>> D=diag([10 10 1]);
>> A=S*D*inv(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.1587e+001
  1.0138e+001
  1.0014e+001
  1.0001e+001
  1.0000e+001
  1.0000e+001
  1.0000e+001
  1.0000e+001
  1.0000e+001
  1.0000e+001
x1 =
 -2.8583e-001
 -9.1466e-001
 -2.8583e-001
niter =
    10

err =
2.2466e+000
2.1028e-001
2.0934e-002
2.0925e-003
2.0924e-004
2.0924e-005
2.0924e-006
2.0924e-007
2.0924e-008
2.0924e-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

function [lambda,v]=power_method(A,x0,maxit)
v=x0;
for 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 11.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.0e+005 *
 -0.8333
 -2.6666

   -0.8333
>> [lambda,v]=power_method(A,x0,100)
lambda =
   10.0000
v =
  1.0e+100 *
   -0.8333
   -2.6667
   -0.8333
>> [lambda,v]=power_method(A,x0,1000)
lambda =
   NaN
v =
   NaN
   NaN
   NaN
>>

La ragione è semplice. Per relativamente piccolo si ha _{}10 _{}e quindi per

10 ^{1} _{}10^{}_{}

da cui

2 10^{}_{}_{2}

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

6.3 Esempio 2

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

\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

0 (1 3) 4 (1 0) 3 ( 1 1)

e quindi il metodo delle potenze applicato ad A, partendo da _{0} può essere sicuramente applicato. D’altra parte dubitiamo converga in quanto _{1} _{2} 1 pur essendo _{1} _{2}. 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
-3.4483e-002
-2.0000e-001
-3.4483e-002
-2.0000e-001
-3.4483e-002
-2.0000e-001
-3.4483e-002
-2.0000e-001
-3.4483e-002
-2.0000e-001
x1 =
3.1623e-001
9.4868e-001
niter =
  16
err =
4.4567e-001
2.4000e+000
4.4567e-001
2.4000e+000
4.4567e-001
2.4000e+000
4.4567e-001
2.4000e+000
4.4567e-001
2.4000e+000
4.4567e-001
2.4000e+000
4.4567e-001
2.4000e+000
4.4567e-001
2.4000e+000
>>

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

6.4 Esempio 3

Per terminare, vediamo il caso della matrice diagonalizzabile (avendo autovalori distinti)

1 2 0 10

in cui il metodo funziona rapidamente, in quanto esiste un solo autovalore di modulo massimo, uguale a 10.

>> A=[1 2; 0 10]
A =
   1     2
   0    10
>> [lambda1, x1, niter, err]=power_basic(A,[1; 3],10^(-8),15)
lambda1 =
9.9779e+000
9.9979e+000
9.9998e+000
1.0000e+001
1.0000e+001
1.0000e+001
1.0000e+001
1.0000e+001
x1 =
2.1693e-001
9.7619e-001
niter =
   8
err =
9.6726e-002
9.7529e-003
9.7610e-004
9.7618e-005
9.7619e-006
9.7619e-007
9.7619e-008
9.7619e-009

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

6.5 Il metodo delle potenze inverse

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

function [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.
% toll: TOLLERANZA.
% nmax: NUMERO MASSIMO DI ITERAZIONI.

%
% OUTPUT:
% lambda : VETTORE DELLE APPROSSIMAZIONI DELL'AUTOVALORE DI MINIMO MODULO.
% x      : AUTOVETTORE RELATIVO ALL'AUTOVALORE DI MINIMO MODULO.
% niter  : NUMERO DI ITERAZIONI.
% err    : VETTORE DEI RESIDUI PESATI RELATIVI A "lambda".
%
% TRATTO DA QUARTERONI-SALERI, "MATEMATICA NUMERICA", p. 184.
%
n=max(size(A)); M=A-mu*eye(n); [L,U,P]=lu(M);
q=z0/norm(z0); q2=q'; err=[]; lambda=[];
res=toll+1; niter=0;
while (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);
    if (costheta > 5e-2)
        res=norm(z-lam*q)/costheta; err=[err; res]; lambda=[lambda; lam];
    else
        disp('\n \t [ATTENZIONE]: AUTOVALORE MULTIPLO'); break;
    end
    x=q;
end

Forniamo ora alcune spiegazioni del codice in invpower.

1. Per risolvere il sistema lineare in 16, si effettua una fattorizzazione della matrice ;

2. All’interno del ciclo while, nella prima riga si calcola _{}, mentre nella successiva un suo versore _{}, e _{}è immagazzinato in lam;

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

6.6 Esempio 1

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

{ \atop }\begin{matrix}15 5 & 7 5 & 1 5 \\ 51 & 25 & 3 \\ 25 5 & 7 5 & 11 5 \end{matrix} { \atop }

{ \atop }1 2 3 2 5 6 7 9 3 { \atop }{ \atop }10 0 0 0 10 0 0 0 1 { \atop }{ \atop }1 2 3 2 5 6 7 9 3 { \atop }^{1}  (21)

Come visto la matrice è 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 del metodo delle potenze consideriamo

0 (1 1 1) (7/6) (1 2 7) 1 (2 5 9) (11/18) (3 6 3)

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

>> 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 11.5]
A =
-15.50000000000000   7.50000000000000   1.50000000000000
-51.00000000000000  25.00000000000000   3.00000000000000
-25.50000000000000   7.50000000000000  11.50000000000000
>> [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 vet-

tore lambda).

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

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

>> s=1-lambda
s =
   0.60983884648007
   0.05762436058732
   0.00573077063146
   0.00057276223344
   0.00005727307685
   0.00000572727622
   0.00000057272730
   0.00000005727272
   0.00000000572727
>> semilogy(1:length(s),s)

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

6.7 Esempio 2: matrice di Poisson

Vogliamo ora applicare il metodo delle potenze e delle potenze inverse per il calcolo di alcuni autovalori di una matrice di Poisson

; ; ; ; ; \begin{matrix} & & 0 & & 0 \\ & & & & 0 \\ 0 & & & & \\ 0 & & & & \\ 0 & 0 & & & \end{matrix} ; ; ; ; ;

con

; ; ; ; ; \begin{matrix}4 & 1 & 0 & & 0 \\ 1 & 4 & 1 & & 0 \\ 0 & 1 & 4 & & \\ 0 & & & & 1 \\ 0 & 0 & & 1 & 4 \end{matrix} ; ; ; ; ;

Salviamo nel file demoautovalori1.m il seguente codice Matlab/Octave:

siz=5;
A = makefish(siz);
eigs=eig(A);
abs_eigs=abs(eigs);
min_eig=min(abs(eigs));
max_eig=max(abs(eigs));
mu=0;
x0=ones(size(A,1),1); nmax=50; toll=10^(-11);
% METODO DELLE POTENZE.
[lambda1, x1, niter, hist_dp]=powerm(A,x0,toll,nmax);
lambdamax=lambda1(length(lambda1));
fprintf('\n \t [MAGGIOR AUTOVALORE (IN MODULO)]');
fprintf('\n \t [CALCOLATO]: %5.15f [ESATTO]: %5.15f',lambdamax,max_eig);
res=norm(A*x1-lambdamax*x1);
fprintf('\n \t [ABS.ERR.]: %2.2e [RES]: %2.2e',
abs(lambdamax-max_eig),res);

% METODO DELLE POTENZE INVERSE.
[lambda, x, niter, hist_ip]=invpower(A,x0,mu,toll,nmax);
lambdamin=lambda(length(lambda));
fprintf('\n \n \t [MINOR AUTOVALORE (IN MODULO)  ]');
fprintf('\n \t [CALCOLATO]: %5.15f [ESATTO]: %5.15f',lambdamin,min_eig);
res=norm(A*x-lambdamin*x);
fprintf('\n \t [ABS.ERR.]: %2.2e [RES]: %2.2e',
abs(lambdamin-min_eig),res);
% METODO DELLE POTENZE INVERSE PER "MU" NON NULLO.
mu=5.9;
[minval,minindex]=min(abs(eigs-mu));
minmu_eig=eigs(minindex);
[lambdamu, x, niter, hist_ipmu]=invpower(A,x0,mu,toll,nmax);
lambdamumin=lambdamu(length(lambdamu));
fprintf('\n \n \t [MINOR AUTOVALORE (IN MODULO)  ]');
fprintf('\n \t [CALCOLATO]: %5.15f [ESATTO]:
%5.15f',lambdamumin,minmu_eig);
res=norm(A*x-lambdamumin*x);
fprintf('\n \t [MU]: %5.5f [ABS.ERR.]: %2.2e [RES]:
%2.2e',mu,abs(lambdamumin-minmu_eig),res);
% METODO DELLE POTENZE INVERSE PER "MU" NON NULLO.
mu2=2.8;
[minval,minindex]=min(abs(eigs-mu2));
minmu2_eig=eigs(minindex);
[lambdamu2, x, niter, hist_ipmu2]=invpower(A,x0,mu2,toll,nmax);

lambdamu2min=lambdamu2(length(lambdamu2));
fprintf('\n \n \t [MINOR AUTOVALORE (IN MODULO)  ]');
fprintf('\n \t [CALCOLATO]: %5.15f [ESATTO]:
%5.15f',lambdamu2min,minmu2_eig);
res=norm(A*x-lambdamu2min*x);
fprintf('\n \t [MU]: %5.5f [ABS.ERR.]: %2.2e [RES]:
%2.2e',mu2,abs(lambdamu2min-minmu2_eig),res);
hold off;
semilogy(1:length(hist_dp),hist_dp,'k-*'); hold on;
semilogy(1:length(hist_ip),hist_ip,'b-o'); hold on;
semilogy(1:length(hist_ipmu),hist_ipmu,'m-d'); hold on;
semilogy(1:length(hist_ipmu2),hist_ipmu2,'r-v');
% legend('DIR','INV','MU1','MU2');

La demo effettua quanto segue.

1. Definita la matrice , i suoi autovalori vengono calcolati con eig, così da poter vedere il comportamento dei vari metodi.

2. Successivamente testiamo il metodo delle potenze per il calcolo dell’autovalore di massimo modulo.

3. In seguito calcoliamo l’autovalore di minimo modulo e quelli più vicini a 5 9 e 2 8.

4. Alla fine plottiamo in scala semilogaritmica, iterazione per iterazione, il residuo (pesato).

5. La parte relativa a

     % legend('DIR','INV','MU1','MU2');

è commentata, in quanto il comando legend non è implementato nelle vecchie versioni di Octave. Per usare la legenda, basta togliere il commento.

Lanciato il programma demoautovalori1, otteniamo

>> demoautovalori1
          [MAGGIOR AUTOVALORE (IN MODULO)]
          [CALCOLATO]: 7.464101615040614 [ESATTO]: 7.464101615137755
          [ABS.ERR.]: 9.71e-011 [RES]: 1.30e-005
          [MINOR AUTOVALORE (IN MODULO)  ]
          [CALCOLATO]: 0.535898384862246 [ESATTO]: 0.535898384862247

Figura 3: Grafico che illustra il residuo per iterazione del metodo delle potenze diretto e nelle versioni di Wielandt.

          [ABS.ERR.]: 1.67e-015 [RES]: 3.46e-012
          [MINOR AUTOVALORE (IN MODULO)  ]
          [CALCOLATO]: 5.732050807568878 [ESATTO]: 6.000000000000000
          [MU]: 5.90000 [ABS.ERR.]: 2.68e-001 [RES]: 1.96e-012
          [MINOR AUTOVALORE (IN MODULO)  ]
          [CALCOLATO]: 2.999999999999130 [ESATTO]: 3.000000000000000
          [MU]: 2.80000 [ABS.ERR.]: 8.71e-013 [RES]: 7.98e-007
>>

Come errore assoluto abbiamo calcolato _{}: _{}dove è l’autovalore esatto mentre _{}è quello calcolato. Quindi abbiamo calcolato la norma 2 del residuo _{2} : dove _{}e l’autovettore calcolato dal metodo delle potenze (diretto o inverso) relativo all’autovalore _{}. Osserviamo che può esserci molta differenza tra _{}e _{2}. Nel metodo delle potenze dirette, e nel calcolo dell’autovalore più vicino a 2 8 si ha infatti _{}_{2}, mentre nel calcolo dell’autovalore più vicino a 5 9 accade il contrario.

6.8 Il metodo QR

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

function [H,Q]=houshess(A)
% REDUCTION OF A MATRIX TO A SIMILAR HESSENBERG ONE.
% SEE QUARTERONI, SACCO, SALERI P. 192.
n=max(size(A)); Q=eye(n); H=A;

for k=1:(n-2)
    [v,beta]=vhouse(H(k+1:n,k)); I=eye(k); N=zeros(k,n-k);
    m=length(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

ove vhouse.m è definito da

function [v,beta]=vhouse(x)
% BUILDING HOUSEHOLDER VECTOR.
% SEE QUARTERONI, SACCO, SALERI P. 197.
n=length(x); x=x/norm(x); s=x(2:n)'*x(2:n); v=[1; x(2:n)];
if (s==0)
    beta=0;
else
    mu=sqrt(x(1)^2+s);
    if (x(1) <= 0)
        v(1)=x(1)-mu;
    else
        v(1)=-s/(x(1)+mu);
    end
    beta=2*v(1)^2/(s+v(1)^2);
    v=v/v(1);
end

quindi QRbasicmethod.m che calcola la matrice triangolare relativa a QR:

function [T,hist]=QRbasicmethod(T_input,maxit)
% QR METHOD FOR A SYMMETRIC TRIDIAGONAL MATRIX "T_input".
T=T_input;
hist=sort(diag(T));
for index=1:maxit
    [Q,R]=qr(T);
    T=R*Q;                     % NEW SIMILAR MATRIX.
    hist=[hist sort(diag(T))]; % IT STORES THE DIAGONAL ELEMENTS
                               % OF THE "index" ITERATION.
end

7 Sugli autovalori della matrice di Poisson

In questa sezione desideriamo calcolare gli autovalori della matrice (simmetrica e a banda) di Poisson che si ottiene dal comando makefish. Più precisamente in demoQR.m scriviamo il seguente codice per il calcolo degli autovalori della matrice di Poisson makefish(2), mostrando come la matrice di partenza viene trasformata in una matrice di Hessenberg da cui si calcola la matrice triangolare superiore prodotta dalle iterazioni di QR:

format short
maxit=100;
siz=2;
A = makefish(siz);
eigs=eig(A);
[H,Q]=houshess(A);
[T,hist]=QRbasicmethod(H,maxit);
A
H
T
eig(A)

Otteniamo come risultato

>> demoQR
A =
     4    -1    -1     0
    -1     4     0    -1
    -1     0     4    -1
     0    -1    -1     4
H =
    4.0000    1.4142    0.0000   -0.0000
    1.4142    4.0000    1.4142   -0.0000
   -0.0000    1.4142    4.0000   -0.0000
         0   -0.0000   -0.0000    4.0000
T =
    6.0000    0.0000   -0.0000   -0.0000
    0.0000    4.0000    0.0000   -0.0000
   -0.0000   -0.0000    4.0000    0.0000
         0   -0.0000    0.0000    2.0000

ans =
2.0000
4.0000
4.0000
6.0000
>>

In cui si vede che gli autovalori, come affermato dalla teoria, sono le componenti diagonali della matrice .

8 Alcuni esempi ed esercizi

1. Lanciare il codice demoQR relativamente alla matrice makefish(3) (e non makefish(2)). A tal proposito modificare la variabile siz.

2. Calcolare con il metodo QR gli autovalori della matrice di Hilbert di ordine 5 (ricordare il comando eig(hilb(5))) e confrontarli con quelli prodotti dal metodo QR. Applicare il metodo delle potenze per il calcolo dell’autovalore di modulo massimo e il metodo di Wielandt per quelle di modulo minimo.

3. Applicare il metodo QR per il calcolo degli autovalori di

{ \atop }0 0 2 1 0 1 0 1 1 { \atop }

E’ a blocchi la matrice generata dal metodo QR? E se cosí non fosse come ovviare al problema? (cf. [6], p. 393). Suggerimento: ricordarsi come calcolare gli autovalori di una matrice di ordine 2 via equazioni di secondo grado.

4. Applicare il metodo QR per il calcolo degli autovalori di

; ; ; ; ; \begin{matrix}3 & 17 & 37 & 18 & 40 \\ 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \end{matrix} ; ; ; ; ;

E’ la matrice generata dal metodo QR a blocchi? E se cosí non fosse come ovviare al problema? (cf. [10], p. 195).

5. Applicare il metodo delle potenze per il calcolo dell’autovalore di massimo modulo di

( ) ; ; ; \begin{matrix} & 2 & 3 & 13 \\ 5 & 11 & 10 & 8 \\ 9 & 7 & 6 & 12 \\ 4 & 14 & 15 & 1 \end{matrix} ; ; ;

Confrontare i risultati ottenuti con quelli di eig(A) per 30 30 plottando il grafico dell’errore tra l’approssimazione della soluzione fornita iterazione per iterazione dal metodo delle potenze con quella di max(abs(eig(A))) (cf. [11], p.133).

6. Si dica se è possibile applicare (sperimentarlo!) il metodo delle potenze per il calcolo dell’autovalore di massimo modulo di

; ; ; \begin{matrix}1/3 & 2/3 & 2 & 3 \\ 1 & 0 & 1 & 2 \\ 0 & 0 & 5/3 & 2/3 \\ 0 & 0 & 1 & 0 \end{matrix} ; ; ;

Dare una spiegazione di tale fenomeno (cf. [11], p. 141, esercizio 6.5).

9 Facoltativo: Metodo QR: alcune considerazioni

Sia una matrice quadrata di ordine e siano da calcolare tutti i suoi autovalori . In precedenza abbiamo visto che é possibile trasformare per similitudine la matrice in una matrice di Hessenberg

; ; ; ; ; ; ; ; \begin{matrix}_{1 1} & _{2 1} & _{3 1} & & _{1} \\ _{1 2} & _{2 2} & _{3 2} & & _{2} \\ 0 & _{2 3} & _{3 3} & & _{3} \\ 0 & 0 & _{3 4} & & _{4} \\ & & & & \\ 0 & 0 & 0 & _{1} & _{} \end{matrix} ; ; ; ; ; ; ; ;

Essendo simile ad , gli autovalori di sono gli stessi di . Si é osservato che una versione di base del metodo QR può incontrare problemi nel convergere qualora esistano 2 autovalori distinti di uguale metodo. Questo é ad esempio il caso della matrice

 ;  ;  ;  ;  ; \begin{matrix}3 & 17 & 37 & 18 & 40 \\ 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \end{matrix}  ;  ;  ;  ;  ;

che ha due valori complessi coniugati

1 27437503036000 1 03039271164378

Affrontiamo il caso della matrice di Hessenberg di ordine 3 (cf. [6], p. 393)

{ \atop }0 0 2 1 0 1 0 1 1 { \atop } Si vede facilmente che alla 25-sima iterazione del metodo QR (in una sua versione di base) corrisponde la matrice

25 { \atop }\begin{matrix}2 000 & 1 069 & 0 926 \\ 0 000 & 0 500 & 0 866 \\ 0 & 0 866 & 0 500 \end{matrix} { \atop }

(abbiamo scritto per semplicità 3 cifre decimali delle componenti di _{25}). Risulta chiaro che quindi il metodo QR non fornisce gli autovalori di , in quanto non tende a una matrice triangolare superiore. Il problema come in precedenza é dovuto alla presenza di 2 autovalori complessi coniugati

0 5 0 8660

Sorge quindi naturale la questione sul come trattare tali casi.

Una seconda osservazione nasce spontanea dal teorema di convergenza del metodo QR. Questa volta consideriamo una sua versione semplificata (cf. [2], p. 625)

Teorema 9.1 Sia una matrice di ordine e supponiamo che i suoi autovalori _{}_{1} soddisfino

1  _{2}  _{}

Allora le iterate _{}del metodo QR convergono a una matrice triangolare superiore che contiene gli autovalori _{}_{1} nelle componenti della diagonale. Se é simmetrica converge a una matrice diagonale. Per quanto riguarda la velocità di convergenza,

max_{}\frac{_{1}}{_{}}(22)

Sono necessarie alcune osservazioni. Rispetto al teorema introdotto in precedenza viene a mancare un’ipotesi sugli autovettori sinistri di che aveva quale effetto il fatto che gli autovalori sulla diagonale di erano ordinati (decrescentemente rispetto al valore delle componenti in modulo). Mentre nel teorema precedente si diceva che le matrici _{}, iterazioni di QR, tendevano a una matrice tridiagonale superiore, qui si parla di _{}ovvero quella matrice triangolare superiore per cui . La cosa comunque basilare é che a guidare la velocità di convergenza é il termine

max_{}\frac{_{1}}{_{}}(23)

Come sfruttare questo fatto per migliorare le proprietà di convergenza del metodo?

10 Facoltativo: Metodo QR con shift

la velocità di convergenza dipende dal rapporto in (23). Se é vicino a 1 la convergenza può essere lenta o perfino non sussistere (si vedano gli esempi della sezione precedente). Per risolvere questi problemi si utilizza una tecnica detta dello shift. Sia un numero (anche complesso!) che approssima meglio degli altri un autovalore e si consideri la successione _{}generata da

\begin{matrix}_{}_{}_{} \\ _{1} _{}_{}_{} \end{matrix}(24)

dove al solito _{}é la matrice identica di ordine , e _{0} : . Osserviamo che essendo _{}matrici unitarie cioé _{}^{}_{}_{}_{}^{}_{}si ha

1 _{}^{}_{}( _{}_{}_{}) _{}^{}_{}_{}_{}_{}

per cui si evince che _{1} ed _{}sono simili. Conseguentemente _{}é simile ad 0 _{}e i suoi autovalori sono _{}_{1}da cui si ottengono facilmente 1. Esistono diverse varianti di questa tecnica.

10.1 Variante 1

In una prima (cf. [3], p. 363) si pone

( _{})_{}1 2  (25)

e si considera il metodo

 ;  ; \begin{matrix}_{}_{}_{}_{} \\ _{1} _{}_{}_{}_{} \\ _{0} \end{matrix}  (26)

Nel caso di matrici hermitiane (cioé ^{}) la convergenza a zero di ( _{})_{1} é del terzo ordine. Una volta che ( _{})_{1} , dove é una tolleranza prefissata, si procede riapplicando lo stesso metodo alla sottomatrice ˆ composta dalle prime 1 righe e colonne di _{}. Poiché ( _{})_{1} 0, e gli autovalori di una matrice diagonale a blocchi sono l’unione degli autovalori di ogni blocco, gli autovalori di sono l’unione di quelli di ˆ con ( _{})_{}.

10.2 Variante 2 (di Wilkinson)

Nel caso in cui _{1} _{}si procede come in precedenza scegliendo _{}l’autovalore della sottomatrice

\begin{matrix}_{1 1} & _{} \\ _{1} & _{} \end{matrix}(27)

che é piú vicino ad ( _{})_{}. Si noti che in questo caso, anche se la matrice ha coefficienti reali, l’utilizzazione dello shift può portare ad una matrice ad elementi complessi, con un aumento del costo computazionale.

11 Online

Sul calcolo degli autovalori di una matrice è possibile trovare online e gratuitamente del software per Matlab.

1. Dal sito

    http://mox.polimi.it/~fausal/matlab.html

si possono scaricare i files eigpower.m e invshift.m che implementano i metodi delle potenze. e qrbasis.m una versione di base del metodo QR.

2. Vario software di algebra lineare è presente al sito

    http://www.math.sc.edu/~meade/math706/MATLAB-files/index.html

in particolare le versioni con shift di QR citate poco sopra.

3. Estremamente interessante la homepage di Y. Saad

    http://www-users.cs.umn.edu/~saad/books.html

in cui si trovano gratuitamente manuali sul calcolo di autovalori e risoluzione di sistemi lineari.

References

[1] E. Andersson e P.A. Ekström, Investigating Google’s PageRank algorithm, http://www.it.uu.se/edu/course/homepage/projektTDB/vt04/projekt5/website/report.pdf .

[2] K. Atkinson, Introduction to Numerical Analysis, Wiley, 1989.

[3] D. Bini, M. Capovani e O. Menchi, Metodi numerici per l’algebra lineare, Zanichelli, 1988.

[4] C. Brezinski e M. Redivo Zaglia, Some numerical analysis problems behind web search, Transgressive Computing 2006.

[5] C. Brezinski e M. Redivo Zaglia, The PageRank vector: properties, computation, approximation and acceleration, SIAM J. Matr. Anal. Appl., Vol.28, N.2, p. 551-575.

[6] V. Comincioli, Analisi Numerica, metodi modelli applicazioni, Mc Graw-Hill, 1990.

[7] S.D. Conte e C. de Boor, Elementary Numerical Analysis, 3rd Edition, Mc Graw-Hill, 1980.

[8] The MathWorks Inc., Numerical Computing with Matlab, http://www.mathworks.com/moler

[9] Netlib, http://www.netlib.org/templates/matlab/

[10] A. Quarteroni, R. Sacco, F. Saleri, Matematica numerica, 2001.

[11] A. Quarteroni e F. Saleri, Introduzione al calcolo scientifico, Springer Verlag, 2006.

[12] A. Suli e D. Mayers, An Introduction to Numerical Analysis, Cambridge University Press, 2003.

[13] Wikipedia (Teoremi di Gerschgorin)

http://it.wikipedia.org/wiki/Teoremi_di_Gerschgorin .

[14] Wikipedia (Inverse iteration)

http://en.wikipedia.org/wiki/Inverse_iteration .

[15] Wikipedia (Metodo delle potenze)

http://it.wikipedia.org/wiki/Metodo_delle_potenze .

[16] Wikipedia (PageRank)

http://it.wikipedia.org/wiki/Page_rank .