AUTOVALORI DI UNA MATRICE^{∗}

A. SOMMARIVA^{†}

Conoscenze richieste. Matrici, vettori. Operazioni con matrici e vettori. Matrici simmetriche. Conoscenze di Matlab/Octave.

Conoscenze ottenute. Teoremi di Gerschgorin. Metodo delle potenze e loro convergenza. Metodo delle potenze inverse e loro convergenza. Metodo QR e loro convergenza.

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

Ax = λx, x ≠ 0  (1.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 n autovalori in ℂ.

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], [34].

Una interessante applicazione è l’algoritmo di PageRank [32], 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], [35]).

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

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

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

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

∗Ultima revisione: 12 febbraio 2010 †Dipartimento di Matematica Pura ed Applicata, Universitá degli Studi di Padova, stanza 419, via Trieste 63,

35121 Padova, Italia (alvise@euler.math.unipd.it). Telefono: +39-049-8271350.

Figura 2.1. Cerchi di Gerschgorin della matrice A definita in (3.2)

Teorema 2.2. (Secondo teorema di Gerschgorin). Se l’unione M_{1} di k cerchi di Gerschgorin è 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 Gerschgorin, dal confronto con la figura 2 abbiamo che 1 autovalore sta nel cerchio K_{3} mentre 2 stanno nell’unione dei cerchi K_{1}, K_{2}.

Teorema 2.3. (Terzo teorema di Gerschgorin). Se la matrice di ordine n è 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 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)}.

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 2.4. Ricordiamo che A è una matrice a coefficienti reali, allora A e A^{T} hanno gli stessi autovalori (cf. [3, p.47]) e quindi applicando i teoremi di Gerschgorin 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

Figura 2.2. Semyon Aranovich Gershgorin (19011933) e Herman Müntz (1884-1956)

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 Gerschgorin ad A^{T} invece che ad A?

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 A una matrice quadrata di ordine n con n autovettori x_{1}, …, x_{n} linearmente indipendenti e autovalori λ_{1}, …, λ_{n} tali che

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

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

1. Una matrice A è diagonalizzabile se e solo se possiede n autovettori linearmente indipendenti.

2. Se tutti gli autovalori di A sono distinti la matrice è diagonalizzabile; l’opposto è ovviamente falso (si pensi alla matrice identica).

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

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

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

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}  (3.3)
y_{k} = Ay_{k−1}, k = 1, 2, …  (3.4)

Teorema 3.1. 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})}  (3.5)

converge a λ_{1}.

Dimostrazione. Per la dimostrazione si confronti [13, 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}

Osserviamo ora che

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

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}. Si osservi che il coefficiente di Rayleigh ρ(·, A) è continuo in ogni x ≠ 0, x ∈ ℝ^{n} in quanto tanto il numeratore quanto il denominatore sono funzioni 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},  (3.7)

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

Nota 3.2. 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 3.3. In virtú di alcune possibili problemi di underflow e underflow si preferisce normalizzare passo passo il vettore y_{k} definito in (3.3). Conseguentemente l’algoritmo del metodo delle potenze risulta

u_{k} = At_{k−1}  (3.8)
t_{k} = \frac{u_{k}}{β_{k}} , β_{k} = ∥u_{k}∥_{2}  (3.9)
l_{k} = ρ(t_{k}, A)  (3.10)

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

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

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

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

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

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

ρ(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}.  (4.4)

da cui è immediato calcolare λ_{n}.

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 per questioni numeriche calcolare λ_{n} via 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 [12, p.181]:

(A − μI) z_{k} = q_{k−1}  (4.5)
q_{k} = z_{k}/∥z_{k}∥_{2}  (4.6)
σ_{k} = (q_{k})^{H}Aq_{k}  (4.7)

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 [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 A una matrice quadrata di ordine n. Utilizzando il metodo di Householder è 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:

1. 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. [3, p.149]) .

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

3. La routine Matlab qr effettua tale fattorizzazione.

Figura 5.1. Helmut Wielandt (1910-2001) e Vera N. Kublanovskaya (1920)

4. 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}.

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}.

Per la convergenza del metodo esistono vari risultati (cf. [7], p. 393, [13, p.180], [11, p.21, Thm. 5.3]). Ricordiamo principalmente [6, p.149]

Teorema 5.1. Se la matrice A ∈ ℝ^{n×n} è regolare e con autovalori tutti distinti in modulo, con

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

Figura 5.2. James Wallace Givens (1910-1993), Karl Hessenberg (1904-1959) e Alston Scott Householder (1904-1993)

allora l’algoritmo QR converge ad una matrice A_{∞} triangolare superiore. Se la matrice è simmetrica, allora

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

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

Alcuni dettagli.

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

2. Per inciso, si può dimostrare che la velocità di convergenza dipende dal rapporto

ρ := max_{1≤i≤n−1} max \frac{|λ_{i+1}|}{|λ_{i}|}.  (5.2)

Il numero di moltiplicazioni necessarie all’algoritmo di Givens per calcolare tale matrice T a partire da A è approssimativamente 10n^{3}/3 mentre per quanto riguarda l’algoritmo di Householder è 5n^{3}/3 [14, p.252]. Il metodo QR applicato ad una matrice A in forma di Hessenberg superiore ha ad ogni passo un costo di 2n^{2} operazioni moltiplicative.

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

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

5.1. Facoltativo: Metodo QR con shift. la velocità di convergenza dipende dal rapporto ρ in (5.2). 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 {A_{k}} generata da

\{ \begin{matrix}A_{k} − μ I = Q_{k} R_{k} \\ A_{k+1} = R_{k} Q_{k} + μ I_{n} \end{matrix}  (5.3)

dove al solito I_{n} é la matrice identica di ordine n, e A_{0} := A. Osserviamo che essendo Q_{k} matrici unitarie cioé Q_{K}^{H} Q_{K} = Q_{K} Q_{K}^{H} = I_{n} si ha

Q_{K} A_{k+1} Q_{K}^{H} = Q_{K} (R_{k} Q_{k} + μ I_{n}) Q_{K}^{H} = Q_{k} R_{K} + μ I_{n} = A_{k}

per cui si evince che A_{k+1} ed A_{k} sono simili. Conseguentemente A_{k} é simile ad A_{0} = A−μ I_{n} e i suoi autovalori sono {λ_{i} −μ}_{i=1,...,n} da cui si ottengono facilmente {λ_{i}}_{i=1,...,n}. Esistono diverse varianti di questa tecnica.

5.1.1. Facoltativo. Variante 1. In una prima variante di QR, (cf. [3], p. 363), detta (a_{k})_{n,n} la componente (n, n) della matrice A_{k}, si pone

μ_{k} = (a_{k})_{n,n} , k = 1, 2, …  (5.4)

e si considera il metodo

\{ \begin{matrix}A_{k} − μ_{k} I = Q_{k} R_{k} \\ A_{k+1} = R_{k} Q_{k} + μ_{k} I_{n} \\ A_{0} = A \end{matrix}  (5.5)

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

5.1.2. Facoltativo. Variante 2 (di Wilkinson). Nel caso in cui |λ_{n−1}| = |λ_{n}| si procede come in precedenza scegliendo μ_{k} l’autovalore della sottomatrice

( \begin{matrix}a_{n−1,n−1} & a_{n,n} \\ a_{n,n−1} & a_{n,n} \end{matrix} )  (5.6)

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

6. Facoltativo: Metodo QR: alcune considerazioni. Sia A una matrice quadrata di ordine n e siano da calcolare tutti i suoi autovalori {λ_{i}}. In precedenza abbiamo visto che é possibile trasformare per similitudine la matrice A in una matrice di Hessenberg T

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

Essendo T simile ad A, gli autovalori di T sono gli stessi di A. 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 i.

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

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

Si vede facilmente che alla 25-sima iterazione del metodo QR (in una sua versione di base) corrisponde la matrice

A_{25} = ( \begin{matrix}2.000 & 1.069 & 0.926 \\ 0.000 & −0.500 & 0.866 \\ 0 & −0.866 & −0.500 \end{matrix} )

(abbiamo scritto per semplicità 3 cifre decimali delle componenti di A_{25}). Risulta chiaro che quindi il metodo QR non fornisce gli autovalori di A, 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 i.

7. Esperimenti numerici in Matlab.

7.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{||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. [12, p.180])

7.2. Esempio 1. Testiamo il codice relativamente al calcolo dell’autovalore 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}  (7.1)

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

7.3. Esempio 2. 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}. 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.

7.4. Esempio 3. Per terminare, vediamo il caso della matrice diagonalizzabile (avendo autovalori distinti)

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

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 toll = 10^{−8}.

7.5. Il metodo delle potenze inverse. Una versione di base invpower del metodo delle potenze inverse [12, 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 4.5, si effettua una fattorizzazione P M = LU della matrice M = A − μI;

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

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

7.6. Esempio 1. 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}  (7.2)

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 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 inverse applicato ad A, e avente quale punto iniziale z_{0} può essere utilizzato per il calcolo dell’autovalore di minimo modulo di A.

>> 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 vettore lambda).

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

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

   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.

7.7. 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 A 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 T 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

8. 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 calcoli il condizionamento della matrice in norma 2 e lo si confronti con cond(hilb(5),2). Eseguire lo stesso esercizio utilizzando il metodo QR.

9. Esercizi facoltativi.

9.1. Facoltativo. Sugli autovalori della matrice di Poisson. In questa sezione desideriamo calcolare gli autovalori della matrice (simmetrica e a banda) di Poisson

A = ( \begin{matrix}B & −I & 0 & … & 0 \\ −I & B & −I & … & 0 \\ 0 & −I & B & … & … \\ 0 & … & … & … & −I \\ 0 & 0 & … & −I & B \end{matrix} )

con

B = ( \begin{matrix}4 & −1 & 0 & … & 0 \\ −1 & 4 & −1 & … & 0 \\ 0 & −1 & 4 & … & … \\ 0 & … & … & … & −1 \\ 0 & 0 & … & −1 & 4 \end{matrix} ) .

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 A viene trasformata in una matrice di Hessenberg H da cui si calcola la matrice T 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 T.

Vogliamo ora applicare il metodo delle potenze e delle potenze inverse per il calcolo di alcuni autovalori di una matrice di Poisson (ottenuta come makefish(5)) i cui blocchi sono matrici 5 × 5. 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);

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

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 A, 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
           [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 e_{a} := |λ−λ_{c}| dove λ è l’autovalore esatto mentre λ_{c} è quello calcolato. Quindi abbiamo calcolato la norma 2 del residuo r_{2} := ∥Ax_{c}−λ_{c}x_{c}∥ dove x_{c} e l’autovettore calcolato dal metodo delle potenze (diretto o inverso) relativo all’autovalore λ ≈ λ_{c}.

Osserviamo che può esserci molta differenza tra e_{a} e r_{2}. Nel metodo delle potenze dirette, e nel calcolo dell’autovalore più vicino a μ = 2.8 si ha infatti e_{a} ≪ r_{2}, mentre nel calcolo dell’autovalore più vicino a μ = 5.9 accade il contrario. In quest’ultimo caso il metodo prima si avvicina all’autovalore 6 e poi diverge verso l’altro autovalore che vale approssimativamente 5.732050807568876e + 00. Quindi il residuo è giusto sia piccolo, ma non si ha la convergenza verso l’autovalore sperato.

9.2. Facoltativo. Calcolo degli zeri di un polinomio, via autovalori della matrice compagna. Si può dimostrare (cf. [29]) che gli zeri del polinomio monico di grado n

p(t) = c_{0} + c_{1}t + ⋯ + c_{n−1}t^{n−1} + t^{n}

sono gli autovalori della matrice compagna (o equivalentemente della sua trasposta visto che eig(A) = eig(A^{T}))

C = [\begin{matrix}0 & 0 & … & 0 & −c_{0} \\ 1 & 0 & … & 0 & −c_{1} \\ 0 & 1 & … & 0 & −c_{2} \\ ⋮ & ⋮ & … & ⋮ & ⋮ \\ 0 & 0 & … & 1 & −c_{n−1} \end{matrix}] .

Sfruttando il comando eig, e osservando che

>> A=[1 2; 3 4]
A =
      1      2

3      4
>> B=[5 6; 7 8]
B =
5      6
7      8
>> C=[A B]
C =
1      2      5      6
3      4      7      8
>> D=[A; B]
D =
1      2
3      4
5      6
7      8
>>

scrivere un codice Matlab che risolva un’equazione algebrica calcolando gli autovalori della matrice compagna associata. Effettuare quindi un test su un polinomio di cui si conoscono gli zeri e valutare la bontà del proprio codice. Di seguito si calcolino gli zeri di modulo massimo e minimo usando il metodo delle potenze e la sua variante inversa.

9.3. Facoltativo: altri esempi ed esercizi. Di seguito citiamo alcuni esercizi ed esempi.

1. Lanciare il codice demoQR relativamente alla matrice definita da makefish(3) (e non come scritto nel file, alla matrice ottenuta con il comando 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

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

E’ a blocchi la matrice generata dal metodo QR? E se cosí non fosse come ovviare al problema? (cf. [7], 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

A = ( \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. [12], p. 195).

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

A(α) = ( \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. [13], p.133).

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

A = ( \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. [13], p. 141, esercizio 6.5).

10. 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.

RIFERIMENTI BIBLIOGRAFICI

[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] C. Brezinski e M. Redivo Zaglia, Méthodes numériques itératives, Ellipses, 2006.

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

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

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

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

[11] P. Olver, Notes on orthogonal bases and the workings of the QR algorithm http://www.math.umn.edu/ olver/aims /qr.pdf

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

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

[14] M. Redivo Zaglia, Calcolo numerico. Metodi e algoritmi, quarta edizione., L. Progetto, 2009.

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

[16] Mac Tutor (Semyon Aranovich Gershgorin)

http://www-history.mcs.st-andrews.ac.uk/Biographies/Gershgorin.html .

[17] Mac Tutor (Alston Scott Householder)

http://www-gap.dcs.st-and.ac.uk/h̃istory/Mathematicians/Householder.html .

[18] Mac Tutor (Helmut Wielandt)

http://www-gap.dcs.st-and.ac.uk/h̃istory/Biographies/Wielandt.html .

[19] NA Digest, John Francis, Co-Inventor of QR http://www.netlib.org/na-digest-html/07/v07n34.html .

[20] Wikipedia (John G.F. Francis)

http://en.wikipedia.org/wiki/John G.F. Francis .

[21] Wikipedia (Semyon Aranovich Gershgorin)

http://en.wikipedia.org/wiki/Semyon Aranovich Gershgorin .

[22] Wikipedia (James Wallace Givens)

http://en.wikipedia.org/wiki/Wallace Givens .

[23] Wikipedia (Karl Hessenberg)

http://en.wikipedia.org/wiki/Karl Hessenberg .

[24] Wikipedia (Alston Scott Householder)

http://it.wikipedia.org/wiki/Alston Scott Householder .

[25] Wikipedia (Vera N. Kublanovskaya)

http://en.wikipedia.org/wiki/Vera N. Kublanovskaya , http://www.pdmi.ras.ru/eng/perso/kublanovskaya.php .

[26] Wikipedia (Herman Müntz)

http://en.wikipedia.org/wiki/Herman Müntz .

[27] Wikipedia (Eigenvalue)

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

[28] Wikipedia (Inverse iteration)

http://en.wikipedia.org/wiki/Inverse iteration .

[29] Wikipedia (Matrice Compagna)

http://it.wikipedia.org/wiki/Matrice compagna .

[30] Wikipedia (Metodo delle potenze)

http://it.wikipedia.org/wiki/Metodo delle potenze .

[31] Wikipedia (PageRank)

http://it.wikipedia.org/wiki/Page rank .

[32] Wikipedia (Rayleigh quotient)

http://en.wikipedia.org/wiki/Rayleigh quotient .

[33] Wikipedia (QR algorithm)

http://en.wikipedia.org/wiki/QR algorithm .

[34] Wikipedia (QR decomposition)

http://en.wikipedia.org/wiki/QR decomposition .

[35] Wikipedia (Teoremi di Gerschgorin) http://it.wikipedia.org/wiki/Teoremi di Gerschgorin .