Maria Grazia Gasparo Metodi numerici per il calcolo di ... · L’algoritmo usato da Google e`...

Maria Grazia Gasparo

Metodi numerici per il calcolo diautovalori e autovettori,valori singolari e vettori singolaridi matrici reali

Indice

1 AUTOVALORI E AUTOVETTORI 11.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Prerequisiti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

1.2.1 Relazione di similitudine e forme canoniche . . . . . . . .41.2.2 Localizzazione degli autovalori . . . . . . . . . . . . . . 81.2.3 Matrici ortogonali: Householder e Givens . . . . . . . . . 9

1.3 Condizionamento del problema . . . . . . . . . . . . . . . . . . . 121.4 Riduzione in forma di Hessenberg . . . . . . . . . . . . . . . . . 171.5 Metodi per matrici non simmetriche . . . . . . . . . . . . . . . . 20

1.5.1 Metodo delle potenze . . . . . . . . . . . . . . . . . . . . 201.5.2 Metodo di Iterazione inversa . . . . . . . . . . . . . . . . 291.5.3 Metodo QR . . . . . . . . . . . . . . . . . . . . . . . . . 31

1.6 Metodi per matrici simmetriche . . . . . . . . . . . . . . . . . . . 441.6.1 Metodo QR per matrici simmetriche tridiagonali . . . . .451.6.2 Metodo Iterazione quoziente di Rayleigh . . . . . . . . . 461.6.3 Metodo Divide–and–conquer . . . . . . . . . . . . . . . . 481.6.4 Metodo di bisezione . . . . . . . . . . . . . . . . . . . . 521.6.5 Metodo di Jacobi . . . . . . . . . . . . . . . . . . . . . . 59

2 DECOMPOSIZIONE AI VALORI SINGOLARI 692.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 692.2 SVD e sue proprieta . . . . . . . . . . . . . . . . . . . . . . . . . 70

2.2.1 Esistenza della SVD . . . . . . . . . . . . . . . . . . . . 702.2.2 SVD e decomposizioni spettrali . . . . . . . . . . . . . . 732.2.3 SVD e approssimazioni di rango basso . . . . . . . . . . 752.2.4 SVD, numero di condizionamento e pseudoinversa . . . . 802.2.5 Rango numerico . . . . . . . . . . . . . . . . . . . . . . 822.2.6 SVD e problema dei minimi quadrati lineari . . . . . . . . 84

2.3 Metodi numerici per il calcolo della SVD . . . . . . . . . . . . . 862.3.1 Caratteristiche generali dei metodi . . . . . . . . . . . . . 862.3.2 Riduzione in forma bidiagonale . . . . . . . . . . . . . . 882.3.3 Metodo QR . . . . . . . . . . . . . . . . . . . . . . . . . 902.3.4 Metodo di Jacobi . . . . . . . . . . . . . . . . . . . . . . 92

Bibliografia 96

1AUTOVALORI E AUTOVETTORI

1.1 Introduzione

La necessita di calcolare alcuni o tutti gli autovalori e/ogli autovettori di una ma-trice nasce in molti problemi pratici, in settori diversi delle scienze applicate. Nellibro di Y. Saad [35] sono ad esempio citati alcuni classici esempi nel campo del-l’ingegneria, della chimica e dell’economia. Nei due esempi seguenti ricordiamoun paio di applicazioni in settori diversi dai precedenti.

Esempio 1.1.1.Il primo esempio riguarda la musica. La fisica degli strumentimusicali e stata ed e oggetto di molti studi. Citiamo a questo proposito il librodi A. H. Benade dal titolo “Fundamental of Musical Acoustics” del 1977, ripub-blicato in una nuova edizione nel 1990 da Dover, il libro “Thephysics of musi-cal instruments” di N.H. Fletcher e T.D. Rossing, pubblicato in 2a edizione nel1998 dalla Springer, il libro di D. Benson dal titolo “Music:a mathematical offer-ing”, comparso nel 1995 e continuamente aggiornato dall’autore, disponibile sulweb all’indirizzohttp://www.maths.abdn.ac.uk/ bensondj/html/maths-music.html,e infine l’articolo divulgativo di V.E. Howle e L.N. Trefethen dal titolo “Eigen-values and musical instruments”, pubblicato nel 2001 (Journal of Computationaland Applied Mathematics, vol. 135, pp. 23-40). In questi studi si creano model-li matematici del “sistema fisico” strumento, solitamente facendo alcune ipotesisemplificative rispetto alla realta, e si studiano le proprieta e le caratteristiche dellostrumento attraverso le proprieta della soluzione del modello. Per moltissimi stru-menti (chitarra, flauto, clarinetto, timpano, ecc.) le frequenze delle oscillazioni,e quindi dei suoni, sono le parti immaginarie degli autovalori di una matrice, eil valore assoluto delle parti reali degli stessi autovalori rappresenta la velocita didecadimento della stessa oscillazione. Introducendo nel modello via via piu in-formazioni sulla realta, ad esempio tenendo conto che una corda di chitarra none di fatto perfettamente flessibile, non e sotto vuoto e quindi le sue vibrazionirisentono dell’aria circostante, e cosı via, si riescono ariprodurre dei grafici degliautovalori che corrispondono sempre di piu al suono effettivamente emesso daglistrumenti. �

Esempio 1.1.2.Un altro campo conosciutissimo in cui gli autovalori e gli au-tovettori giocano un ruolo fondamentale e il campo dei motori di ricerca sul web,

2 Capitolo 1

in particolare Google. Quando inseriamo una o piu parole chiave per una ricerca,Google seleziona le pagine contenenti le parole nella sua banca dati e rispondealla richiesta presentando queste pagine in ordine di importanza, o autorevolezza.L’autorevolezza di una paginap dipende da quante pagine la citano (ovvero con-tengono un link ap) e da quanto queste pagine sono a loro volta autorevoli. Questoconcetto viene utilizzato per la selezione e l’ordinamentodelle pagine medianteun algoritmo PageRank, inventato alla fine degli anni ’90 daifondatori di Google,Sergey Brin e Lawrence Page. L’algoritmo usato da Google e sicuramente moltocomplesso ed e un segreto industriale. Nonostante cio, esiste una notevole mole diarticoli che studiano possibili algoritmi PageRank. Senzaentrare nei dettagli, percui rimandiamo alla letteratura (per esempio al delizioso articolo di K. Bryan eT. Leise, “ The $25,000,000,000 Eigenvector: the linear algebra behind Google”,SIAM Review 48, 2006,pp. 569-581), accenniamo alle basi delpiu semplice al-goritmo PageRank che si possa concepire, per scoprire che questo coinvolge gliautovalori di un’opportuna matrice. SiaP l’insieme delle pagine presenti nel webattinenti alle parole chiave inserite, en la sua cardinalita, e facciamo l’ipotesi sem-plificativa che tutte le pagine diP contengano almeno un link ad un’altra paginadi P stesso. In questa ipotesi possiamo definire l’autorevolezza xk dellak-esimapagina mediante la formula

xk =∑

j∈Lk

xj

nj, (1.1)

doveLk e l’insieme degli indici delle pagine che contengono un link alla paginak e nj e il numero totale di link a pagine contenute inP che laj-esima pagi-na contiene. Supponiamo in questa formula di non contare eventuali link di unapagina a se stessa. Con la (1.1) affermiamo il principio che una pagina che fariferimento allak conta di piu se contiene pochi links, fra cui quello che ci inte-ressa, perche e una pagina “selettiva”, mentre conta menose contiene molti links.La formula (1.1) non e utilizzabile direttamente per calcolare l’autorevolezza diuna singola pagina. Basti pensare al caso non infrequente che due pagine, peresempio lar-esima e las-esima, si citino a vicenda; in questa situazione la for-mula diventa inutilizzabile perche per calcolarexr devo conoscerexs e viceversa.D’altra parte la (1.1) puo essere utilizzata per calcolarel’autorevolezza di tutte len pagine coinvolte. Scrivendo infatti simultaneamente len relazioni e indicandoconx il vettore(x1, . . . , xn)T , si ottiene un sistema della forma

Ax = x, (1.2)

doveA e una matricen× n, detta “matrice dei links”, con elementi dati da

aij =

{

1nj

se la paginaj rimanda alla paginai0 altrimenti.

Possiamo osservare che la (1.2) ammette soluzione diversa da quella banale, cheovviamente non ci interessa, se e solo seλ = 1 e autovalore diA, nel qual casoxe un autovettore corrispondente. A questo proposito osserviamo che nelle ipotesi

AUTOVALORI E AUTOVETTORI 3

fatte la matriceA e stocastica per colonne, ovvero ha tutti elementi non negativie la somma degli elementi su ciascuna colonna e uguale a 1. Daquesto segue cheAT e = e, dove abbiamo indicato cone il vettore con tutte le componenti uguali a1. Qesta ultima relazione dice cheλ = 1 e autovalore diAT e quindi anche diA. Ilfatto cheA ammetta l’autovalore uguale a 1 non basta ai fini pratici; occorre ancheche questo autovalore sia semplice, in modo che l’autospazio corrispondente abbiadimensione 1 e non ci siano difficolta nello scegliere il vettore delle autorevolezze.Si puo far vedere cheλ = 1 e sicuramente semplice se len pagine sono tuttecollegate fra loro, ovvero se da una qualunque di esse e possibile raggiungerneun’altra qualunque di link in link, e se non ci sono pagine prive di links. Purtroppoquesta situazione e abbastanza improbabile in una banca dati come quella del webfatta da bilioni di pagine. Fortunatamente si puo modificare opportunamente lamatrice dei link in modo cheλ = 1 sia effettivamente un suo autovalore semplice.�

Come e ben noto, gli autovalori di una matriceA sono in generale numericomplessi e, se la parte immaginaria di un autovalore e diversa da zero, gli au-tovettori corrispondenti sono anch’essi complessi. Soltanto nel caso di matricireali e simmetriche si ha la sicurezza che autovalori e autovettori sono reali. Perquesto motivo, in molti testi di analisi numerica (come per esempio [5] o [40])i metodi per il calcolo di autovalori e/o autovettori sono descritti per matrici inC

n×n, cioe a elementi nel campo dei numeri complessi. Noi preferiamo limitar-ci al caso di matrici reali, anche se questo in alcuni casi complica la trattazione,per due ordini di motivi. Primo, perche in molte applicazioni le matrici che siincontrano sono reali. Secondo, perche non tutti gli studenti a cui ci rivolgiamohanno molta familiarita con i numeri complessi. Nel seguito di questo capitolo,con il termine “numero complesso” intenderemo un numero appartenente aC\R,ossia un numero complesso con parte immaginaria diversa da zero. Altrimenti, ilnumero sara detto “reale”.

Sappiamo che gli autovalori di una matriceA ∈ Rn×n sono le radici del poli-

nomio caratteristicop(λ) = det(A − λI) e i metodi per calcolarli sono pertantometodi iterativi, a meno che non sian = 2. Di solito i metodi sono organizzatiin due fasi. In un primo tempo si trasforma la matriceA in una matrice simileBdotata di una struttura che faciliti, o renda molto piu economico, il calcolo effetti-vo degli autovalori/autovettori. La struttura scelta e quella di Hessenberg nel casogenerale e tridiagonale nel caso simmetrico, a cui si arrivamediante un numerofinito di trasformazioni diA effettuate tramite matrici ortogonali1. Ci occupere-mo di questa trasformazione nel paragrafo 1.4. La seconda fase e quella iterativa,nella quale si applica un metodo iterativo appropriato allastruttura della matrice(simmetrica o no) e allo scopo che ci prefiggiamo (calcolare un autovalore, alcuniautovalori, tutti gli autovalori, con autovettori corrispondenti o meno). Nel para-

1Fa eccezione il metodo di Jacobi, di cui ci occuperemo nel paragrafo 1.6.5, che si applica amatrici simmetriche senza trasformarle preventivamente in forma tridiagonale.

4 Capitolo 1

grafo 1.5 introdurremo i metodi per matrici generali (possiamo pensarle semprein forma di Hessenberg): partiremo dal metodo delle potenzeche sotto opportuneipotesi permette di calcolare l’autovalore di modulo massimo e un autovettore adesso associato; passeremo poi al metodo di iterazione inversa che viene utilizza-to per calcolare un autovettore corrispondente ad un autovalore approssimato conqualche altro metodo; approderemo poi al metodo principe, noto come metodoQR, con il quale si calcolano tutti gli autovalori. Ovviamente questi metodi pos-sono essere utilizzati anche per matrici simmetriche, ma esistono metodi specificiche sfruttano la simmetria diA, di cui parleremo nel paragrafo 1.6. Oltre al meto-do di Jacobi, vedremo altri due metodi, entrambi applicabili a matrici tridiagonalie simmetriche: il primo di questi va sotto nomi diversi (metodo di Sturm, meto-do di Householder, metodo di bisezione) ed e essenzialmente un’applicazione delmetodo di bisezione alla ricerca di uno, alcuni, o tutti gli zeri del polinomio carat-teristico; l’altro, noto come metodo divide–and–conquer,viene usato per calcolaretutti gli autovalori e i corrispondenti autovettori.

Per quanto riguarda il software, la funzioneeig di Matlab 2 per il calcolodi autovalori e/o autovettori usa il metodo QR per qualsiasimatrice. Invece, lalibreria di sottoprogrammi open source per l’algebra lineare LAPACK (disponi-bile in linguaggio FORTRAN e in C) prevede anche i metodi di bisezione edivide–and–conquer per matrici simmetriche (vedi [1] e www.netlib.org).

Notazioni. Nel seguito useremo i seguenti simboli:

• I per la matrice identita; quando dovesse risultare necessario specificarne ladimensione useremo un pedice (ad esempio,I3 sara la matrice identita di di-mensione 3).

• ei per l’i-esimo vettore della base canonica diRn qualunque sian; quando ne-

cessario, specificheremo la dimensione dello spazio con un apice (ad esempio,em1 sara il primo vettore della base canonica diR

m).• Λ(A) per lo spettro di una matriceA.• xH per il trasposto coniugato dix, siax un vettore o una matrice.

1.2 Prerequisiti

1.2.1 Relazione di similitudine e forme canoniche

Ricordiamo che due matriciA e B sono simili se esiste una terza matriceSinvertibile tale che

A = SBS−1.

In questo casoA eB hanno gli stessi autovalori, mentre gli autovettori diA sonoottenuti da quelli diB previa moltiplicazione perS. La maggior parte dei meto-di numerici per il calcolo di autovalori/autovettori si basano su trasformazioni

2Qui e nel seguito faremo riferimento alla versione Matlab 7.8.0.


di similitudine, nel senso che traggono le loro basi teoriche dall’esistenza di al-cuneforma canonicheche, tramite un’opportuna relazione di similitudine, fanno“emergere” gli autovalori. Elenchiamo alcune fra le principali forme canoniche acui faremo riferimento nel seguito del capitolo.

Diagonalizzazione.Una matriceA si dicediagonalizzabilese e simile a una ma-trice diagonale, ossia se esistonoX,Λ ∈ C

n×n, conX invertibile eΛ diagonale,tali che

A = XΛX−1.

In questo casoΛ contiene sulla diagonale gli autovalori diA e le colonne diXcostituiscono un sistema completo di autovettori linearmente indipendenti.Ricordiamo cheA e diagonalizzabile se e soltanto se e non difettosa3 e chel’insieme delle matrici diagonalizzabili e denso inR

n×n.

Diagonalizzazione unitaria, o ortogonale.Una matriceA si diceunitariamentediagonalizzabilese esistonoQ,Λ conQ unitaria (ortogonale nel caso reale)4 eΛdiagonale, tali che

A = QΛQH . (1.3)

Ricordiamo che l’insieme delle matrici unitariamente diagonalizzabili coincidecon l’insieme delle matrici normali (ossia quelle matrici per cui vale l’uguaglianzaAAH = AHA), di cui fanno parte ad esempio le matrici reali simmetriche.

Forma di Schur. Un risultato fondamentale dell’algebra lineare dice che qualsiasimatriceA ∈ R

n×n puo essere fattorizzata nella forma

A = QTQH , (1.4)

conQ e T in Cn×n, Q unitaria eT triangolare superiore. La (1.4) e nota come

decomposizione di Schur, o forma di Schur, diA. Da questa segue cheA e similealla matriceT ed ha pertanto gli stessi autovalori, che peraltro compaiono sulladiagonale diT stessa. SeA non e simmetrica e ammette autovalori complessi,alloraT eQ sono complesse. Ma seA e simmetrica, allora le due matrici sonoreali e in particolareT deve essere simmetrica e pertanto diagonale; allora lafattorizzazione di Schur (1.4) coincide con la diagonalizzazione ortogonale diA.

Qui di seguito dimostriamo il teorema di Schur, che stabilisce l’esistenza del-la fattorizzazione (1.4) per qualsiasi matrice quadrata, anche in campo complesso.

3Si dice che una matrice edifettosa quando qualche autovalore ha molteplicita algebricamaggiore della molteplicita geometrica.

4Ricordiamo che una matriceQ ∈ Cn×n e unitaria seQHQ = QQH = I. Una matrice reale

unitaria e detta ortogonale.

6 Capitolo 1

Teorema 1.2.1.DataA ∈ Cn×n, esistono due matriciQ e T in C

n×n, conQunitaria eT triangolare superiore, tali cheA = QTQH .

Dimostrazione. La dimostrazione procede per induzione sun. Pern = 1, la tesie vera conT = A eQ = (1). Pern > 1, siaλ un autovalore diA, x un autovet-tore normalizzato corrispondente aλ, e{y2, . . . , yn} una base ortonormale per ilcomplemento ortogonale diS = span{x}. Allora, la matriceQ1 = (x y2 · · · yn)e unitaria e tale cheQH

1 x = e1. Consideriamo ora la matriceB = QH1 AQ1; la

prima colonna diB e data da

Be1 = QH1 AQ1e1 = QH

1 Ax = λQH1 x = λe1,

e pertantoB ha la forma a blocchi seguente:

B =

(

λ cH

0 A1

)

.

Per ipotesi induttiva la matriceA1 ∈ C(n−1)×(n−1) ammette una forma di Schur

A1 = Q2T2QH2 , conQ2 unitaria eT2 triangolare superiore. Allora

B =

(

λ cH

0 Q2T2QH2

)

,

da cui segue

A = Q1BQH1 = Q1

(

λ cH

0 Q2T2QH2

)

QH1 .

Ponendo

Q3 =

(

1 0H

0 Q2

)

si ottiene

A = Q1Q3

(

λ cHQ2

0 T2

)

QH3 Q

H1 = Q1Q3TQ

H3 Q

H1 ,

ovveroA = QTQH , doveQ = Q1Q3 e unitaria in quanto prodotto di matriciunitarie eT e triangolare superiore .

Forma di Schur reale. E possibile dimostrare [21] che ogni matrice reale am-mette unaforma di Schur realein cui T eQ sono inR

n×n e T e triangolare ablocchi, ovvero ha la forma

T =

T11 T12 · · · T1m

0 T22 · · · T2m...

.... . .

...0 0 · · · Tmm

,


dove ogni blocco diagonaleTii e 1 × 1 reale o2 × 2 con autovalori complessiconiugati. E evidente che lo spettro diT , che coincide con quello diA, contienetutti gli elementiTii di dimensione 1 e gli autovalori deiTii di dimensione 2. Ilmetodo QR, citato nell’introduzione come il principale metodo per il calcolo degliautovalori, e di fatto finalizzato al calcolo della forma diSchur reale.

Forma di Schur e sottospazi invarianti.La forma di Schur e strettamente legataal concetto di sottospazio invariante, che andiamo a definire. SiaS un sottospaziodi R

n; S e invariante per la matriceA se l’operatoreA portaS in se, ovvero

x ∈ S ⇒ Ax ∈ S.

L’esempio piu semplice di sottospazio invariante eS = span{x}, dovex e unautovettore diA; questo sottospazio ha ovviamente dimensione 1. In generale,seS ha dimensionep e {x1, . . . , xp} e una sua base, chiamiamoX la matricein Rn×p che hax1, . . . , xp come colonne; allora possiamo far vedere cheS esottospazio invariante diA se e solo se esiste una matriceB ∈ Rp×p tale che

AX = XB. (1.5)

Infatti, da xi ∈ S segueAxi ∈ S, ovvero esiste un vettoreyi ∈ Rp tale che

Axi = Xyi; da questo segue immediatamente la (1.5) ponendoB = (y1 · · · yp).Viceversa, siax un generico elemento diS, che puo essere scritto comeXy perqualchey ∈ R

p; alloraAx = AXy = XBy = Xz e pertantoAx ∈ S.

Da (1.5) si deduce una relazione molto interessante fraΛ(A) eΛ(B). Infatti,siaλ ∈ Λ(B) ey ∈ Rp un autovettore corrispondente. Allora, daBy = λy segueAXy = XBy = λXy, ovveroλ ∈ Λ(A) eXy e un autovettore corrispondente.In altri termini,Λ(B) ⊂ Λ(A). Si dice cheS e il sottospazio invariante associatoa questo sottoinsieme diΛ(A).

Per vedere il legame fra la forma di Schur (1.4) e i sottospaziinvarianti diA,fissiamo un interop ∈ {1, . . . , n} e poniamo

Q = (Q1 Q2) e T =

(

T11 T12

0 T22

)

,

conQ1 ∈ Rn×p eT11 ∈ Rp×p. Allora la (1.4) diventa

A(Q1 Q2) = (Q1 Q2)

(

T11 T12

0 T22

)

,

da cui segue in particolareAQ1 = Q1T11;

questa relazione dice che le colonne diQ1 costituiscono una base di un sot-tospazio invariante diA di dimensionep e cheΛ(T11) ⊂ Λ(A). Inoltre, daAQ2 = Q1T12 + Q2T22 segue che seA e simmetrica (nel qual casoT12 = 0)anche le colonne diQ2 generano un sottospazio invariante.

8 Capitolo 1

1.2.2 Localizzazione degli autovalori

Poniamoci la seguente domanda. Data una matrice realeA non diagonale, in qualezona del piano complesso (o dell’asse reale nel caso cheA sia simmetrica) sonodisposti gli autovalori? La risposta piu immediata si ha ricordando che il raggiospettrale di una matrice e l’estremo inferiore delle normenaturali della matricestessa, da cui si evince che

Λ(A) ⊂ {z ∈ C : |z| ≤ ||A||},

dove|| · || e una qualsiasi norma naturale, ad esempio norma–1 o norma–∞.

Una risposta meno grossolana puo essere ottenuta dalle seguenti considera-zioni. Consideriamo un autovaloreλ e un autovettorex ad esso associato, e siar l’indice della componente di massimo modulo dix. Allora dall’uguaglianza∑n

j=1 arjxj = λxr segue

λ− arr =

n∑

j=1,j 6=r

arjxj

xr

e quindi

|λ− arr| ≤n∑

j=1,j 6=r

|arj |.

Pertantoλ ∈ {z ∈ C : |z − arr| ≤∑n

j=1,j 6=i |arj |}. Non sappiamo a prioriquale e l’indicer, ma possiamo comunque concludere che

λ ∈ R = ∪ni=1Ri,

dove

Ri = {z ∈ C : |z − aii| ≤n∑

j=1,j 6=i

|aij|}.

Facendo le stesse considerazioni con la matriceAT , si deduce anche che

λ ∈ C = ∪nj=1Cj,

dove

Cj = {z ∈ C : |z − ajj| ≤n∑

i=1,i6=j

|aij |}.

Un teorema che comprende e precisa le considerazioni ora fatte e il seguente teo-rema, dovuto a S.A. Gerschgorin (1901-1933), la cui dimostrazione puo esseretrovata ad esempio in [19].


Teorema 1.2.2.DataA ∈ Rn×n, consideriamo l’insiemeR = ∪n

i=1Ri. Allora,ad ogni componente connessa massimale diR appartengono tanti autovalori diA quanti sono i cerchi compresi nella componente stessa, contando autovalori ecerchi con le loro molteplicita. Analogo risultato vale per l’insiemeC.

I cerchiR1, . . . , Rn sono detti “cerchi di Gerschgorin per righe” eC1, . . . , Cn

sono detti “cerchi di Gerschgorin per colonne”.

1.2.3 Matrici ortogonali: Householder e Givens

In dimensione 2 ci sono soltanto due tipi di matrici ortogonali. Il primo tipo ecostituito dalle matrici di riflessione, matrici simmetriche della forma

(

−c ss c

)

, (1.6)

dovec = cosθ es = sinθ per un dato angoloθ. La generalizzazione di queste ma-trici a dimensionen generica porta alle matrici di Householder, che sono matricisimmetriche della forma

U = I − 2uuT

‖u‖2, (1.7)

doveu e un qualsiasi vettore diverso da zero inRn. E facile vedere che (1.6) e

una matrice di Householder corrispondente au = (cosθ2 , sinθ

2)T .Le matrici di Householder sono detteriflettori elementaria causa del loro

significato geometrico. Siau fissato; consideriamo un vettorex e scomponiamo-lo nelle sue componenti lungo la direzione diu e quella perpendicolare, ovveroscriviamox = αu + w, conα scalare ew ∈ R

n tale chewTu = 0. ApplicandoU adx otteniamo

Ux = α(I − uuT

β)u+ w − 1

βuuTw = αu− α

βuuTu+ w = −αu+ w.

Il vettoreUx ha quindi la stessa componente dix lungo la direzione perpendico-lare au, mentre la componente lungo la direzione diu risulta uguale ma di segnoopposto: in praticaU ha riflessox rispetto alla direzione diw.

Le matrici di Householder vengono in generale usate per annullare “un pez-zo” di un vettore datox ∈ R

n. Ad esempio, si puo facilmente verificare chescegliendo

u = x± ‖x‖2e1, (1.8)

si ottieneUx = ∓‖x‖2e1,

ovvero si annullano tutte le componenti dix esclusa la prima. Un problema nu-merico importante in questa situazione e la scelta del segno nella definizione diu.Tenendo conto cheu1 = x1 ± ‖x‖2, una scelta comune e

u1 = x1 + sign(x1)‖x‖2, (1.9)

10 Capitolo 1

che permette di evitare errori di cancellazione numerica quando|x1| ≃ ‖x‖2 (insostanza, quandox e vicino ad essere un multiplo die1). Questa non e l’uni-ca scelta possibile; ad esempio, in [21] si propone un’altrascelta utile quando sivuole cheUx sia un multiplo positivo die1; in questo caso si deve usare neces-sariamente il segno “−” nella (1.8) e per evitare eventuali errori di cancellazionesi puo utilizzare la seguente formula:

u1 =−(x2

2 + . . . + x2n)

x1 + ‖x‖2. (1.10)

La costruzione precedente puo essere generalizzata nel modo seguente. Fis-sato un interok ∈ {0, . . . , n− 1}, poniamo

u = (0, . . . , 0, xk+1±s, xk+2, . . . , xn)T , con s =√

x2k+1 + · · · + x2

n. (1.11)

Si verifica facilmente che la matrice di HouseholderU = I − 2 uuT

‖u‖2 agisce suxnel modo seguente:

Ux = (x1, . . . , xk,∓s, 0, . . . , 0)T . (1.12)

In altri termini,U lascia intatte le primek componenti dix, ne cambia la(k+ 1)-esima e annulla le ultimen − k − 1. Data la forma diu, e immediato vedere cheU ha una forma a blocchi:

U =

( Ik 0

0 U

)

, (1.13)

doveU e la matrice di Houseolder di dimensionen − k che agisce sulle ultimen− k componenti dix annullandole tutte meno che la prima.

Una proprieta che rende le matrici di Householder molto comode da usarenella pratica e la seguente: per calcolare il prodotto fra la matriceU e un vettorex non occorre costruire esplicitamenteU , ma basta conoscereu; infatti, ponendoβ = uT u

2 , si ha

Ux =

(

I − 1

βuuT

)

x = x− 1

βuuTx = x− uTx

βu. (1.14)

Tornando in dimensione 2, il secondo tipo di matrici ortogonali e rappresen-tato dalle matrici di rotazione, matrici non simmetriche della forma

U =

(

c s−s c

)

, (1.15)

dove ancorac = cosθ e s = sinθ per un dato angoloθ. Il nome “matrici dirotazione” deriva dal fatto che il vettorey = Ux risulta ruotato rispetto ax di


un angoloθ in senso orario. Per verificare questa asserzione, esprimiamo x incoordinate polari, ad esempiox = ρ(cosφ, sinφ)T ; calcolando il vettorey si ot-tieney = ρ(cos(φ− θ), sin(φ− θ))T . E facile vedere che la matriceUT provocaanch’essa una rotazione di un angoloθ, ma in senso antiorario.

Spesso avremo la necessita di utilizzare una rotazione perannullare una com-ponente di un dato vettore; piu in particolare, vorremo scegliere θ in modo chesia

(

c s−s c

)T (ab

)

=

(

∗0

)

, (1.16)

per dati valori dia e b, ovveroas + bc = asinθ + bcosθ = 0. E interessanteosservare che questa equazione puo essere risolta senza calcolare esplicitamenteθ. A titolo di esempio riportiamo l’algoritmo seguente, tratto da [21], dove sicalcolano direttamentec es in funzione diτ = tgθ o di τ = cotgθ, scegliendo frale due possibilita quella che evita l’insorgere di possibili fenomeni di overflow.

Algoritmo 1.2.1. Algoritmo per risolvere il problema (1.16)Dati: a, bRisultati: c, sSeb = 0, allora: c = 1 es = 0

altrimenti:se|b| > |a|, allora:

τ = −ab,s = 1√

1+τ2, c = sτ

altrimenti:τ = − b

a,c = 1√

1+τ2, s = cτ

Fine sceltaFine scelta

La generalizzazione delle matrici di rotazione (1.15) al cason× n porta allematrici di Givens, che sono ottenute dalla matrice identit`a cambiando soltanto glielementi di indici(p, q), (q, p), (p, p) e (q, q), conp, q ∈ {1, · · · , n} diversi fraloro fissati. D’ora in poi supporremo semprep < q. Una volta fissatip eq, si fissaanche un angoloθ, e si definisce la matrice di GivensG = G(p, q, θ) mediante lerelazioni seguenti:

gij =

{

1 se i = j0 se i 6= j

,per i, j 6= p, q

gpp = gqq = cgpq = −gqp = s,

dove, come prima,c = cosθ e s = sinθ. La struttura di queste matrici puo essereampiamente sfruttata quando sono coinvolti prodotti fra matrici. Ad esempio, sevogliamo sovrascrivere a una matriceB ∈ R

n×n il prodottoGTB possiamo usareil seguente semplice algoritmo, in cui si tiene conto del fatto cheGT modificasoltanto le righe di indicip e q di B:

12 Capitolo 1

Algoritmo 1.2.2. Algoritmo per il calcolo diGTBPerj = 1, . . . , n

t1 = bpj

t2 = bqj

bpj = ct1 − st2bqj = st1 + ct2

Fine ciclo

Analogamente, l’algoritmo seguente sovrascrive aB il prodottoBG:

Algoritmo 1.2.3. Algoritmo per il calcolo diBGPeri = 1, . . . , n

t1 = bipt2 = biqbip = ct1 − st2biq = st1 + ct2

Fine ciclo

1.3 Condizionamento del problema

Lo studio del condizionamento dello spettro di una matrice puo essere affrontatoin modo efficace esaminando un autovalore per volta. Cominciamo dal caso di unautovalore semplice. Sia pertantoλ un autovalore semplice diA, eA + δA unamatrice perturbata, con‖δA‖ piccola. Sia inoltreλ+ δλ l’autovalore diA+ δAcorrispondente aλ, ovvero l’autovalore della matrice perturbata piu vicinoa λ.Ci chiediamo se e possibile stimare|δλ| in funzione di‖δA‖ 5. Una risposta aquesta domanda e fornita dal seguente teorema nel quale si suppone senza perditadi generalita cheδA abbia la formaδA = ǫE, doveǫ e uno scalare (piccolo) eE ∈ R

n×n e una matrice fissata.

Teorema 1.3.1.Sia λ un autovalore semplice diA, x e y i corrispondenti au-tovettori normalizzati destro e sinistro6 e E ∈ R

n×n una matrice qualunque.Allora esiste una funzioneλ(ǫ) analitica per|ǫ| sufficientemente piccolo tale cheλ(0) = λ, λ(ǫ) e autovalore semplice diA+ ǫE e

|λ(ǫ) − λ| ≤ ‖E‖|yHx|ǫ+ O(ǫ2). (1.17)

5Nel campo degli autovalori lo studio del condizionamento viene fatto in generale in senso as-soluto, a differenza di altri settori della matematica, neiquali invece si considerano le perturbazionirelative. Un trattamento esaustivo dell’argomento puo essere trovato in [38]. Un accenno alla teoriadel condizionamento in senso relativo puo invece essere trovato in [8].

6Un autovettore sinistroy perA e un autovettore (destro) perAT ; pertanto soddisfa la relazioneyHA = λyH .


Dimostrazione. Premettiamo che la (1.17) ha senso perche per un autovaloresemplice il prodottoyHx e sempre diverso da zero (cfr. [26, Lemma 6.3.10]).L’esistenza diλ(ǫ) analitica tale cheλ(0) = λ e λ(ǫ) e autovalore semplice diA + ǫE discende da argomenti di teoria delle funzioni di variabilecomplessa,come si puo trovare ad esempio in [39, Teorema 6.8.8]. Indichiamo ora conx(ǫ)un autovettore relativo aλ(ǫ). Derivando rispetto aǫ la relazione

(A+ ǫE)x(ǫ) = λ(ǫ)x(ǫ)

si ottieneEx(ǫ) + (A+ ǫE)x(ǫ) = λ(ǫ)x(ǫ) + λ(ǫ)x(ǫ)

e, perǫ = 0,

Ex+Ax(0) = λ(0)x+ λx(0).

Moltiplicando a sinistra peryH otteniamo

λ(0) =yHEx

yHx

e possiamo pertanto scrivere i primi termini dello sviluppodi λ(ǫ)

λ(ǫ) = λ+yHEx

yHxǫ+ O(ǫ2),

da cui ricaviamo la (1.17).

In pratica il teorema ci dice che un autovalore semplice puoessere mal con-dizionato (in senso assoluto) se

s =1

|yHx|

e grande. In altri termini,s funge da numero di condizionamento.Il teorema appena dimostrato vale per perturbazioni sufficientemente piccole.

In realta, risultati analoghi (nello spirito) possono essere dimostrati senza questarestrizione, a patto di trovare maggiorazioni meno stringenti di quella espressadalla (1.17). A titolo di esempio riportiamo un teorema di Bauer-Fike, la cuidimostrazione puo essere trovata in [8].

Teorema 1.3.2.Supponiamo cheA abbia tutti autovalori sempliciλ1, λ2, . . . , λn

e sianoxi e yi autovettori normalizzati destro e sinistro relativi aλi. Sia inoltreδA una perturbazione introdotta sulla matriceA. Allora lo spettro diA + δA einteramente contenuto nell’unioneB1 ∪ B2 ∪ · · · ∪ Bn, doveBi e il cerchio nelpiano complesso di centroλi e raggion ‖δA‖

|yHi xi|

.

14 Capitolo 1

I risultati precedenti sono assolutamente generali. Come corollario di questirisultati si deduce chegli autovalori semplici di matrici simmetriche sono ben con-dizionatidal momento che la simmetria diA implica l’uguaglianza di autovettoridestri e sinistri e quindiyH

i xi=1 per ognii. Non solo, ma l’ipotesi di simmetria etalmente potente da implicare il buon condizionamento degli autovalori indipen-dentemente dalla loro molteplicita. In questa direzione va il seguente teorema diWeyl, la cui dimostrazione si puo trovare in [8].

Teorema 1.3.3.SianoA e A = A + E due matrici simmetriche. Siano inoltreλ1, λ2, . . ., λn gli autovalori diA (non necessariamente distinti) eλ1, λ2, . . ., λn

quelli di A. Allora|λi − λi| ≤ ‖E‖2

per i = 1, 2, . . . , n.

Torniamo al caso generale, nel quale un autovalore semplicepuo essere (molto)mal condizionato se|yHx| e (molto) piccolo, ovvero sex e y sono (molto) viciniad essere ortogonali. Vediamo qualche esempio in proposito.

Esempio 1.3.1.Dato un numero reale positivoα, consideriamo la matrice

A =

(

0 1α 0

)

,

i cui autovalori sonoλ1 =√α e λ2 = −√

α. Esaminiamo il condizionamen-to di λ1. E facile vedere che gli autovettori destro e sinistro normalizzati sonorispettivamente

x =1√

1 + α

(

1√α

)

e y =1√

1 + α

( √α

1

)

,

cosı cheyHx = 2√

α1+α

. Quandoα e molto piccolo,yHx si avvicina a zero el’autovalore risulta mal condizionato. Analogo ragionamento si puo fare perλ2.�

Nell’esempio appena fatto, i due autovalori tendono a coincidere quandoαtende a zero. Questo potrebbe far pensare che il cattivo condizionamento possaessere necessariamente legato alla vicinanza di due autovalori. In realta non ecosı, come il successivo esempio, dovuto a J.H. Wilkinson,mostra.

Esempio 1.3.2.SiaA ∈ Rn×n definita da:

aij =

{

i per i = jn per i = j − 10 altrimenti .

Gli autovalori sono1, . . . , n e sono pertanto tutti semplici e “ben separati” fraloro. Pern = 10 abbiamo calcolato con la funzioneeig di Matlab gli autovettori


destri e sinistri e i prodotti scalaris1 = yH1 x1, s2 = yH

2 x2, . . . , sn = yHn xn otte-

nendo i seguenti valori:

s1 = · · · = s6 = 0s7 = 6.7763e− 21s8 = 5.4210e− 20s9 = 8.1315e− 20s10 = 1.0842e− 19.

I primi sei valori in teoria dovrebbero essere diversi da zero, ma evidentementesono talmente piccoli da risultare zero in macchina; i successivi quattro valorisono al di sotto della precisione di macchinaǫm ≃ 2.22e−16. In conclusione,questi valori testimoniano un condizionamento estremamente cattivo di tutti gliautovalori. E infatti, se perturbiamo la matrice introducendo nella posizione(n, 1)il valore

√ǫm gli autovalori diventano:

0.99996, 2.0004, 2.9985, 4.0035, 4.99486.0052, 6.9965, 8.0015, 8.9996, 10.000

con perturbazioni rispetto agli autovalori originali che arrivano all’ordine di10−3.Se ripetiamo lo stesso esperimento pern = 20, la maggior parte degli autovaloridella matrice perturbata sono addirittura complessi:

2.0264e+ 17.3576e− 119.497 ± i1.133617.937 ± i2.522915.803 ± i3.720413.259 ± i4.79797.7408 ± i4.51865.1965 ± i3.72043.0633 ± i2.52291.5026 ± i1.1336

�

L’analisi del condizionamento di autovalori multipli per matrici non simme-triche e abbastanza complessa e richiede strumenti di algebra lineare che esulanodai prerequisiti medi dei lettori a cui ci rivolgiamo. Chi einteressato puo co-munque trovare tutti i risultati nel libro [45] di J.H. Wilkinson. Qui osserviamosoltanto che in alcuni casi non e possibile dimostrare un risultato analogo a (1.17),ma si puo dimostrare che|λ(ǫ) − λ| dipende da una potenza frazionaria diǫ. Peresempio, consideriamo la matrice

A =

(

1 10 1

)

16 Capitolo 1

che ha un autovalore doppioλ = 1, conx = (1 0)T e y = (0 1)T , da cui segueyTx = 0; se perturbiamo l’elementoa21 portandolo da zero aǫ, gli autovaloridiventano1 ±√

ǫ.

Torniamo al teorema di Weyl, che sancisce il buon condizionamento degliautovalori di matrici simmetriche.E importante osservare che questo ottimo con-dizionamento e inteso in senso assoluto e quindi non implica che gli autovaloridella matrice perturbata siano affetti tutti dallo stesso errore relativo. Per precisareulteriormente questo concetto, consideriamo il caso particolare in cuiA = fl(A).In questo caso possiamo affermare che la matrice di perturbazioneE soddisfa ladisuguaglianza

‖E‖2 ≤ nǫm‖A‖2. (1.18)

Infatti sappiamo cheeij = aij(1 + ǫij) con |ǫij| < ǫm e pertanto|eij | ≤ ǫm|aij |per ogni coppia di indicii e j. Allora ‖E‖1 ≤ ǫm‖A‖1. Sfruttando l’equivalenzafra norma–1 e norma–2

1√n‖X‖1 ≤ ‖X‖2 ≤ √

n‖X‖1

si ottiene la relazione (1.18).Sia oraλmax l’autovalore diA con il massimo valore assoluto; allora sappi-

amo che‖A‖ = |λmax| e la (1.18) diventa pertanto

‖E‖2 ≤ nǫm|λmax|;

dal teorema di Weyl segue che per ognii ∈ {1, . . . , n} vale

|λi − λi||λi|

≤ n|λmax||λi|

ǫm ≡ kiǫm.

L’accuratezza relativa con cuiλi approssimaλi dipende quindi dalla costanteki:se|λi| e piccolo rispetto a|λmax|, alloraki e grande eλi ha poche cifre significa-tive. Al limite, senkiǫm e maggiore di 1,λi non ha neanche una cifra buona. Inaltre parole, perche ci possa essere un minimo di accuratezza occorre che

|λi| ≥ n|λmax|ǫm. (1.19)

Possiamo quindi concludere che in aritmetica floating point, indipendentementedalla bonta dell’algoritmo utilizzato, gli autovalori piccoli rispetto a quello dimassimo valore assoluto saranno in generale approssimati con scarsa accuratezzarelativa.


1.4 Riduzione in forma di Hessenberg

Una matriceH ∈ Rn×n e detta di Hessenberg superiore se ha la forma

H =

∗ ∗ ∗ ∗ · · · ∗ ∗∗ ∗ ∗ ∗ · · · ∗ ∗

∗ ∗ ∗ · · · ∗ ∗∗ ∗ · · · ∗ ∗

. . . . . .. . . . . .

∗ ∗

.

ossia sehij = 0 peri = 3, . . . , n ej = 1, . . . , i−2. La matrice e dettairriducibilese gli elementi della codiagonale inferiore sono tutti diversi da zero. Spesso nelseguito avremo a che fare con matrici di Hessenberg e supporremo che siano ir-riducibili. Questa ipotesi non lede la generalita del problema del calcolo degliautovalori. Infatti, seH avesse un elemento codiagonale uguale a zero, avrebbela forma triangolare a blocchi

H =

(

H1 Z0 H2

)

Poiche il determinante di una matrice triangolare a blocchi e il prodotto dei de-terminanti dei blocchi diagonali (fatto facilmente dimostrabile per induzione), siotterrebbeΛ(H) = Λ(H1) ∪ Λ(H2). In altri termini, il problema del calcolodegli autovalori di una matrice di Hessenberg non irriducibile si disaccoppia inpiu problemi per matrici di Hessenberg irriducibili.

Qualsiasi matriceA ∈ Rn×n e ortogonalmente simile a una matrice di Hes-

senberg, nel senso che esiste una matriceQ ortogonale tale che

QTAQ = H, (1.20)

conH di Hessenberg superiore. Questa affermazione viene dimostrata costrutti-vamente, facendo ricorso ad esempio alle matrici di Householder. Piu in partico-lare, possiamo prendere matrici della forma (1.13) conk = 1, . . . , n− 2, di voltain volta scegliendo in modo opportuno il vettore x di cui annullare gli elementi dal(k + 2)-esimo all’n-esimo. Al primo passo si prende comex la prima colonnadi A, si costruisce una matriceU1 con la regola (1.13) conk = 1 e si calcolaU1AU1. PoicheU1 e simmetrica e ortogonale, si ottiene una matrice simile adAche ne conserva pertanto gli autovalori. Tenendo conto della (1.13), si vede chel’applicazione diU1 a sinistra lascia invariata la prima riga diA. La successivaapplicazione diU1 a destra lascia invariata la prima colonna diU1A e ne modifica

18 Capitolo 1

le colonne successive. In conclusione, abbiamo ottenuto una matrice della forma

U1AU1 =

∗ ∗ ∗ ∗ · · · ∗ ∗∗ ∗ ∗ ∗ · · · ∗ ∗0 ∗ ∗ ∗ · · · ∗ ∗0 ∗ ∗ · · · ∗ ∗0

. . . . . .

0. . . . . .

0 ∗ ∗

.

Ora consideriamok = 2. Prendiamo come vettorex la seconda colonna diU1AU1, costruiamo la matriceU2 con la (1.13) e calcoliamoU2U1AU1U2. Sipuo facilmente vedere che la moltiplicazione a sinistra per U2 lascia invariate leprime due righe diU1AU1 e non modifica gli zeri gia introdotti nella prima colon-na. La successiva moltiplicazione a destra perU2 non tocca le prime due colonnedi U2U1AU1 e modifica le successive. In conclusione, si ottiene una matrice dellaforma:

U2U1AU1U2 =

∗ ∗ ∗ ∗ · · · ∗ ∗∗ ∗ ∗ ∗ · · · ∗ ∗0 ∗ ∗ ∗ · · · ∗ ∗0 0 ∗ ∗ · · · ∗ ∗0 0

. . . . ..

0 0. .. .. .

0 0 ∗ ∗

.

Successivamente, costruiamoU3 in modo da annullare gli ultimin − 4 elementidella terza colonna diU2U1AU1U2, e cosı via. Dopon − 2 passaggi avremoottenuto

H = Un−2 · · ·U2U1AU1U2 · · ·Un−2 (1.21)

ovvero, ponendoQ = U1U2 · · ·Un−2, la (1.20).Descriviamo ora mediante l’algoritmo 1.4.1 il procedimento appena tratteggia-

to. Nell’algoritmo non abbiamo definito i vettoriu ∈ Rn secondo la (1.11): di vol-

ta in volta, il vettorev rappresenta invece il sottovettore diu con il quale definirela matriceU della (1.13). Abbiamo inoltre deciso di normalizzare il vettore v,cosı cheU = I − 2vvT . L’algoritmo restituisce le matriciH e Q; quest’ulti-ma potrebbe servire qualora si volessero calcolare gli autovettori di A a partireda quelli diH. In alternativa alla costruzione esplicita diQ potremmo resti-tuire come risultati tutti i vettoriv, diciamov1, . . . , vn−2, che contengono tuttele informazioni necessarie a ricostruireQ (cfr. [21, algoritmo 7.7.1]). Inoltre, sipotrebbero definire i vettoriv con la regola (1.10) invece che con la (1.9). Comespecificato in [21], questo algoritmo ha un costoO(n3).


Algoritmo 1.4.1. Algoritmo per la riduzione in forma di Hessenberg superiorediA ∈ R

n×n.Dati: ARisultati: H di Hessenberg superiore eQ ortogonale tali cheQTAQ = H1.H = A2.Q = I3. Perk = 1, . . . , n − 2

1. v = H(k + 1 : n, k)2. v(1) = v(1) + sign(v(1))‖v‖2

3. v = v/‖v‖2

4. H(k + 1 : n, k : n) = H(k + 1 : n, k : n) − 2v(vTH(k + 1 : n, k : n))5. H(1 : n, k + 1 : n) = H(1 : n, k + 1 : n) − 2(H(1 : n, k + 1 : n)v)vT

6. Q(1 : n, k + 1 : n) = Q(1 : n, k + 1 : n) − 2(Q(1 : n, k + 1 : n)v)vT

Fine ciclo4. Fine

Esempio 1.4.1.Se applichiamo l’algoritmo precedente alla matrice

A =

1 1 1 18 4 2 127 9 3 164 16 4 1

,

operando in Matlab si ottiene

H =

1.0000 −1.4159 −0.9931 −0.0947−69.9214 5.3710 16.8125 7.8182

0.0000 1.9479 2.8401 050360 −0.0000 −0.4546 −0.2111

,

doveh31 = 3.5527e−15 eh42 = −8.8818e−16, e

Q =

1.0000 0 0 00 −0.1144 −0.7858 −0.60780 −0.3861 −0.5285 0.75600 −0.9153 0.3212 −0.2430

,

con‖QTAQ−H‖ ≃ 8.35e− 15. Se invece applichiamo l’algoritmo alla matricesimmetrica

A =

1 1 1 11 2 3 41 3 6 101 4 10 20

,

si ottiene

H =

1.0000 −1.7321 −0.0000 −0.0000−1.7321 20.6667 10.2740 −0.0000

0.0000 10.2740 7.1754 −0.36460.0000 0.0000 −0.3646 0.1579

20 Capitolo 1

con elementi “nulli” dell’ordine di10−15 o 10−16, e

Q =

1.0000 0 0 00 −0.5774 0.6556 0.48670 −0.5774 0.0937 −0.81110 −0.5774 −0.7493 0.3244

,

con‖QTAQ−B‖ ≃ 4.78e−15. Si noti cheH e di fatto una matrice tridiagonalee simmetrica. �

Il fatto che la seconda matriceH nell’esempio sia tridiagonale non e un caso.Infatti, dalla (1.20) segue che, seA e simmetrica, ancheH lo e. Ma una matricedi Hessenberg simmetrica e necessariamente tridiagonale(e simmetrica). Cosı,sotto l’ipotesi di simmetria diA, l’algoritmo da una matrice della forma:

α1 β1 0 0 · · · 0β1 α2 β2 0 · · · 00 β2 α3 β3 · · · 0

. . . . . . . . .. . . . . . βn−1

0 0 0 0 βn−1 αn

. (1.22)

Si potrebbe a prima vista pensare che, con un procedimento analogo a quelloappena visto, sia possibile diagonalizzare una matriceA in un numero finito dipassaggi. Sappiamo che questo pensiero non e corretto: se infatti fosse vero,avremmo trovato un metodo diretto per trovare le radici di unpolinomio di gradon qualsiasi! Dove sta il trucco? La risposta e semplice. Per diagonalizzareAdovremmo partire conk = 0 nella (1.13), prendendo comex la prima colonnadi A. In questo modo la moltiplicazioneU0A annullerebbe gli elementiai,1 peri = 2, . . . , n, e questo ci sta bene. Ma la successiva moltiplicazione a destraperU0 rimodificherebbe gli elementi della prima colonna diU0A, probabilmentereintroducendo numeri diversi da zero al posto degli zeri appena introdotti. Inaltre parole la matriceU0AU0 non avrebbe gli zeri nella prima colonna comedesiderato!

1.5 Metodi per matrici non simmetriche

1.5.1 Metodo delle potenze

Data una matriceA qualsiasi, preferibilmente in forma di Hessenberg, il metodoprincipe oggi usato per calcolare tutti gli autovalori e ilmetodoQR. Per capirebene questo metodo occorre partire da un altro metodo, piu antico, il metodo dellepotenze, la cui idea risale al lavoro [30] di C.L. Muntz del 1913. Questo metodosi basa su alcune proprieta delle potenze di matrici che andiamo a descrivere.


Supponiamo cheA sia diagonalizzabile, con autovaloriλ1, . . . , λn tali che

|λ1| > |λ2| ≥ |λ3| ≥ · · · ≥ |λn| (1.23)

e un sistema di autovettori indipendentix1, x2, . . . , xn. L’autovaloreλ1, che esicuramente reale, viene detto autovaloredominante. Sia orav0 6= 0 un qualsiasivettore inR

n, che pertanto puo essere espresso come combinazione lineare degliautovettori, diciamo

v0 = α1x1 + α2x2 + · · · + αnxn.

A partire dav0 definiamo la successione{vk} mediante la regola

vk+1 = Avk = Akv0.

Sfruttando il fatto cheAxi = λixi peri = 1, . . . , n, si vede chevk =∑n

i=1 αiλki xi

per ognik. Seα1 6= 0, possiamo anche scrivere

vk = α1λk1

[

x1 +

n∑

i=2

(

αi

α1

)(

λi

λ1

)k

xi

]

, (1.24)

da cui deduciamo, in virtu della (1.23), che quandok tende all’infinito i vettorivk tendono a disporsi nella direzione dix1. Osserviamo subito che l’insiemedei vettori per i qualiα1 e uguale a zero ha misura nulla inRn, e pertanto laprobabilita di sceglierev0 in tale insieme e praticamente nulla, se ci si affida a unascelta casuale. D’altra parte, se anche dovessimo scegliere v0 in questo insieme,cio non rappresenterebbe un problema in precisione finita,perche gli errori diarrotondamento farebbero spuntar fuori una componente nonnulla lungox1, senon inv0, sicuramente piu avanti. Una volta tanto, gli errori di arrotondamentosono di aiuto invece che di intralcio! Supponiamo quindi nelseguito, senza perditadi generalita,α1 > 0. Sgombrato il campo da questo problema, occupiamoci dialtre due questioni. La prima questione e la seguente. La formula (1.24) mostrache

limk→∞‖vk‖ = α1‖x1‖limk→∞|λ1|k =

{

0 se|λ1| < 1∞ se|λ1| > 1α1‖x1‖ se|λ1| = 1

Poiche in generale|λ1| 6= 1, la formula (1.24) porterebbe in pratica ad underflowo overflow. Poiche siamo interessati alla direzione dell’autovalorex1, possiamorisolvere questo problema calcolando una successione di vettori normalizzati

vk =Avk−1

‖Avk−1‖. (1.25)

22 Capitolo 1

In alcuni (vecchi) libri di testo si consiglia di usare la norma ∞ perche il suocalcolo non richiede operazioni aritmetiche, ma in generale e preferibile usare lanorma euclidea per motivi che saranno chiari fra breve. Posto

yk =n∑

i=2

(

αi

α1

)(

λi

λ1

)k

xi,

si ha

vk =α1λ

k1

|α1λk1 |

(x1 + yk)

‖x1 + yk‖= µk

x1 + yk

‖x1 + yk‖,

doveµk e uguale a 1 per ognik seλ1 > 0, mentre oscilla fra+1 e−1 seλ1 < 0.Pertanto la successione{vk} converge a x1

‖x1‖ nel primo caso e oscilla tra duepunti di accumulazione,x1

‖x1‖ e − x1

‖x1‖ , nel secondo caso. D’altra parte, seλ1 enegativo, entrambi i punti di accumulazione risultano autovettori corrispondentia λ1, e quindi il calcolo della successione ci serve comunque a individuare ladirezione dell’autovettore corrispondente all’autovalore dominante. Abusando deltermine, diremo talvolta nel seguito che la successione{vk} converge a± x1

‖x1‖ .La seconda questione rimasta in sospeso e la seguente. Si possono estendere

le idee descritte finora per riuscire a calcolare una successione di scalari conver-gente aλ1? La risposta e positiva, e lo strumento matematico da utilizzare e ilquoziente di Rayleigh. Ricordiamo che, data una matriceA e un vettorex 6= 0, ilquoziente di Rayleigh associato e

r(x) =xTAx

xTx.

Sex e un autovettore diA corrispondente ad un autovaloreλ, e immediato vederecher(x) = λ. Allora, per ognik possiamo calcolare

σk = r(vk) = (vk)TAvk, (1.26)

dove l’ultima uguaglianza segue dal fatto chevk ha norma euclidea unitaria (ed ec-co perche in fase di normalizzazione conviene usare la norma euclidea). Essendor(x) un funzionale continuo, si ha

limk→∞σk = limk→∞µ2k

(x1 + yk)T

‖x1 + yk‖A

x1 + yk

‖x1 + yk‖= r(x1) = λ1.

In altri termini, anche se la successione dei vettorivk non dovesse convergere insenso stretto all’autovettore, la successione deiσk converge sempre aλ1.

Il metodo delle potenzeconsiste nel calcolo delle successioni{vk} e {σk}mediante le formule (1.25) e (1.26) fino a convergenza della seconda. L’algoritmo1.5.1 realizza il metodo, proseguendo i calcoli fino al soddisfacimento di un crite-rio di arresto o all’esaurimento dellekmax iterazioni consentite. Come criterio diarresto possiamo usare un criterio sul residuo

‖Avk − σkvk‖ ≤ η|σk| (1.27)


oppure un criterio (misto) sulla vicinanza di due successive stime dell’autovalore

|σk − σk−1| ≤ η(|σk−1| + 1). (1.28)

Nel caso si usi questo secondo criterio occorre tenere in memoria due stime con-secutive dell’autovalore, aggiornandole ad ogni iterazione (vedi le variabiliσv eσ nell’algoritmo).

Algoritmo 1.5.1. Algoritmo delle potenze per calcolare l’autovalore dominantedi una matriceA ∈ R

n×n, e un autovettore corrispondente.Dati: n,A, η, kmax, vRisultati: ind indicatore del risultato

ind = 0: il criterio di arresto e stato soddisfattoind = 1: kmax iterazioni sono risultate insufficienti

k numero di iterazioni effettuateσ ultima approssimazione dell’autovalore dominante calcolatav ultima approssimazione dell’autovettore corrispondente

1. v = v‖v‖

2. w = Av3. σv = vTw2. Perk = 1, . . . , kmax

1. v = w‖w‖

2. w = Av3. σ = vTw4. Se il criterio di arresto e soddisfatto, allora:ind = 0, Fermati

altrimenti:σv = σFine ciclo suk

3. ind = 14. Fine

Per quanto riguarda la velocita di convergenza, possiamo dire cheσk tendea λ1 almeno con la velocita con cuiyk tende a zero; si tratta quindi di una con-vergenza lineare con costante pari al rapporto|λ2|

|λ1| : piu l’autovaloreλ1 e domi-nante rispetto agli altri autovalori, e piu veloce e la convergenza; viceversa, laconvergenza puo essere molto lenta quando|λ1| e |λ2| sono molto vicini.

Esempio 1.5.1.Dato un vettorey = (y1, . . . , yn)T , consideriamo la matrice diVandermondeA definita da

aij = yn−ji .

Pern = 5 e y = (1, 2, 3, 4, 5)T , gli autovalori diA sono

45.9604,−28.9437, 6.7954,−0.8496 e 0.0375.

Il rapporto∣

∣

∣

λ2

λ1

∣

∣

∣ e circa 0.63. Partendo dal vettorev = (1, 1, 1, 1, 1)T e usando

la tolleranzaη = 10−12, l’algoritmo 1.5.1 impiega 58 iterazioni per soddisfare

24 Capitolo 1

il criterio di arresto (1.27) fornendo in uscita un’approssimazione dell’autovaloredominante affetta da un errore relativo pari a3.5 × 10−12, mentre ne richiede 63per soddisfare il criterio (1.28).

Se si costruisce invece la matrice partendo dal vettorey = (.1, .2, .3, .4, .5)T ,gli autovalori sono

1.7551,−0.2991, 0.0467,−0.0048 e 0.0002,

con un rapporto∣

∣

∣

λ2

λ1

∣

∣

∣circa uguale a 0.1704. In questo caso, partendo dal solito

vettorev, l’algoritmo impiega 15 e 16 iterazioni, rispettivamente,per soddisfare idue criteri di arresto, con un errore relativo sull’autovalore calcolato di4.2×10−13

nel primo caso e7.2 × 10−14 nel secondo. �

Le considerazioni fatte sulla velocita di convergenza delmetodo delle poten-ze prescindono dalla struttura diA. Di fatto, un’analisi piu attenta mostra chela velocita e piu elevata nel caso cheA sia simmetrica. Infatti, consideriamo ilquozienter(x) e calcoliamone il gradiente. Con semplici manipolazioni alge-briche si ottiene

∇r(x) =(xTx)(A+AT )x− 2(xTAx)x

(xTx)2,

da cui, seA e simmetrica,

∇r(x) =2(Ax− r(x)x)

xTx.

Allora, ∇r(x) = 0 se e soltanto sex e un autovettore e, di conseguenza,r(x)il corrispondente autovalore. Supponiamo ora per comodit`a ‖x1‖ = 1. Comeabbiamo gia osservato, seλ1 > 0, la successione{vk} converge ax1 e possiamoscrivere

σk − λ1 = r(vk) − r(x1)= ∇r(x1)

T (vk − x1) + O(‖vk − x1‖2)= O(‖vk − x1‖2).

Seλ1 < 0 si possono ripetere gli stessi passaggi separatamente per l’insiemedegli indiciK tali che la sottosuccessione{vk}k∈K converge ax1 e per l’insiemecomplementare degli indiciK per cui la sottosuccessione{vk}k∈K converge a−x1. In ogni caso, poiche

‖vk − (±x1)‖ = O(

∣

∣

∣

∣

λ2

λ1

∣

∣

∣

∣

k)

.

si arriva alla conclusione che seA e simmetrica

σk − λ1 = O(

∣

∣

∣

∣

λ2

λ1

∣

∣

∣

∣

2k)

.


Abbiamo descritto e studiato il metodo delle potenze nell’ipotesi che la ma-triceA sia diagonalizzabile e ammetta un autovalore dominante. Sipuo facilmentevedere che il metodo funziona anche nel caso che l’autovalore dominante abbiamolteplicitaq maggiore di 1, ovvero quando

λ1 = λ2 = . . . = λq

e|λ1| > |λq+1| ≥ · · · ≥ |λn|,

purche la molteplicita geometrica sia anch’essaq (ipotesi sicuramente soddisfattaseA e diagonalizzabile). In questo caso la successione{vk} converge ad unvettore appartenente al sottospazio generato dax1, . . . , xq. La convergenza non einvece assicurata quando gli autovaloriλ1, λ2, . . . , λq hanno lo stesso modulo masono distinti (ad esempio,q = 2, λ1 eλ2 complessi coniugati).

Esempio 1.5.2.Consideriamo un simpatico esempio che nasce dalla geometriaelementare, trattato nel recente articolo [17], a cui rimandiamo per tutti i dettagliche per questioni di spazio qui trascureremo. Datin punti casuali nel piano carte-siano, consideriamo il poligono che ha questi punti come vertici. Indichiamo conx0 e y0 i vettori che raccolgono le ascisse e le ordinate dei punti rispettivamente,e conP0 = P (x0, y0) il poligono. A partire daP0 costruiamo un nuovo poligonoP1 = P (x1, y1) i cui vertici sono i punti di mezzo dei lati diP0, indicando conx1 e y1 i vettori che raccolgono le ascisse e le ordinate dei nuovi vertici. E facilevedere chex1 = Mx0 ey1 = My0, doveM e la matrice

M =1

2

1 1 0 0 · · · 00 1 1 0 · · · 00 0 1 1 · · · 0

. . . . . . . . .. . . . . . 1

1 0 0 0 0 1

.

A partire daP1 costruiamo ora con le stesse modalita di prima un nuovo poligonoP2 = P (x2, y2), dovex2 = Mx1 e y2 = My1. E cosı via, generiamo unasuccessione di poligoni

Pk = P (xk, yk),

dovexk = Mxk−1 e yk = Myk−1.

La domanda e: dove ci porta questa costruzione? Cosa e il limite dei poligoni perk tendente all’infinito, ammesso che questo limite esista? L’intuizione geometricaci dice che i poligoni dovrebbero tendere a un punto. L’analisi numerica confermaquesta intuizione. Consideriamo a questo scopo la successione {xk}. Questasuccessione e ottenuta applicando il metodo delle potenze(senza normalizzazioni)

26 Capitolo 1

alla matriceM a partire dal vettorex0. Allora, seguendo i passaggi fatti perarrivare alla relazione (1.24), possiamo affermare che

xk = α1λk1

[

z1 +n∑

i=2

(

αi

α1

)(

λi

λ1

)k

zi

]

, (1.29)

doveλ1, . . . , λn sono gli autovalori diM ordinati in senso decrescente di modulo,z1, . . . , zn sono i corrispondenti autovettori eα1, α2, . . . , αn sono le componentidi x0 lungo gli autovettori stessi; per poter ricavare la (1.29) dobbiamo supporreα1 = zT

1 x0 6= 0. Analogamente, per la successione{yk} possiamo scrivere

yk = β1λk1

[

z1 +

n∑

i=2

(

βi

β1

)(

λi

λ1

)k

zi

]

, (1.30)

doveβ1, β2, . . . , βn sono le componenti diy0 lungo gli autovettori e supponiamoβ1 = zT

1 y0 6= 0.Per poter dire qualcosa di preciso sul comportamento delle due successioni

occorre studiare gli autovalori e autovettori diM . A questo scopo osserviamocheM = 1

2(I + S), doveS = (en e1 e2 . . . en−1) e una permutazione dellamatrice identita. Gli autovaloriω1, . . . , ωn di S sono noti e dati dalle radicin-esime dell’unita

ωj = cos

(

2πj

n

)

+ isin

(

2πj

n

)

, j = 1, . . . , n

e i corrispondenti autovettori normalizzativ1, . . . , vn sono

vj =

√

1

n

1ωj

ω2j...

ωn−1j

, j = 1, . . . , n.

Gli autovalori diM sono pertanto12(1 + ωj), per j = 1, . . . , n, e gli autovet-tori restanov1, . . . , vn. Possiamo ordinare gli autovalori diM considerandoli nelseguente modo:

λ1 = 12(1 + ω0) = 1

λ2 = 12(1 + ω1)

λ3 = 12(1 + ωn−1) = λ2

λ4 = 12(1 + ω2)

λ5 = 12(1 + ωn−2) = λ4

. . .

,

da cui segue (supponiamon dispari per semplicita) che

1 = |λ1| > |λ2| = |λ3| > · · · > |λn−1| = |λn|.


Visto il riordino fatto per gli autovalori, gli autovettoridiM corrispondenti diven-tano:

z1 = v1 = e√n

z2 = v1z3 = vn−1

z4 = v2z5 = vn−2

. . .

,

dovee indica il vettore con tutte le componenti uguali a 1. Ora che abbiamo questeinformazioni, possiamo tornare alla (1.29) e concludere che, essendoλ1 = 1, lasuccessione{xk} converge al vettoreα1z1, senza alcun bisogno di normalizzarei vettori ad ogni iterazione. Inoltre, essendoα1 = zT

1 x0, dall’espressione diz1deduciamo che

xk → eTx0

ne = ξe

doveξ = 1n

∑ni=1 x0,i ex0,i e l’i-esima componente dix0 (ossia l’ascissa dell’i-

esimo vertice del poligono inizialeP0).Analogamente, per la successione{yk} si deduce che

yk → eT y0

ne = ηe

con η = 1n

∑ni=1 y0,i. Dal punto di vista geometrico questo risultato dice che i

vertici dei poligoniPk tendono a convergere tutti a un unico punto di coordinate(ξ, η); questo punto altro non e che il centroide dei vertici diP0.

La velocita con la quale i vertici dei poligoni collassano verso il centroidedipende dal rapporto fra i moduli diλ1 eλ2, dati da

|λ1| = 1 e |λ2| =

√

1 + cos(2πn

)√

2.

Evidentemente, il modulo diλ2 tende a 1 al crescere din, ed e abbastanza vicino a1 gia per valori din non troppo grandi, come indicato nello specchietto seguente:Da questo deduciamo che la convergenza sara tanto piu lenta quanto piu numerosi

n 5 7 · · · 17 19 · · · 25 · · · 55|λ2| 0.8090 0.9010 0.9830 0.9864 0.9921 0.9984

sono i punti scelti inizialmente.

Tutti i ragionamenti ora fatti valgono seα1 e β1 sono entrambi diversi dazero. Nell’articolo [17] ci si chiede cosa succede seα1 = β1 = 0, ovvero se ilcentroide dei vertici diP0 e nell’origine degli assi. In questo caso, dallo sviluppoxk = α2λ

k2z2 +α3λ

k3z3 + · · ·+αnλ

knzn e dal fatto che tutti gli autovalori in esso

28 Capitolo 1

presenti hanno modulo minore di 1, deduciamo che la successione{xk} convergeal vettore nullo. La stessa cosa vale per{yk}.

Proviamo ad applicare allora il metodo delle potenze vero e proprio, nel qualei vettori vengono normalizzati ad ogni iterazione. Cosa succede in questo caso?La teoria del metodo delle potenze dice che, essendoci due autovalori dominanticon lo stesso modulo ma diversi fra loro il metodo non converge. In effetti, in [17]si dimostra analiticamente che la successione dei poligoniPk = P (xk, yk) ha due“poligoni di accumulazione”P e ¯P ; al primo tendono i poligoni di indice pari eall’altro tendono i poligoni di indice dispari. Si dimostrainoltre che i vertici dientrambi questi poligoni appartengono a un’ellisse di cui si ricava esplicitamentel’equazione.

−1 −0.5 0 0.5−1

−0.5

0

0.5poligono 0

−1 −0.5 0 0.5−1

−0.5

0

0.5poligono 9

−0.5 0 0.5−0.5

0

0.5poligono 30

−0.5 0 0.5−0.5

0

0.5poligono 31

Figura 1.1 Esempio 1.5.2: dal poligono iniziale P0 ai due “poligoni diaccumulazione” P e ¯P

Per verificare questo risultato, abbiamo fatto diverse prove con diversinscegliendo dei vettorix0 ey0 casuali e traslati in modo che il centroide dei verticidi P0 = P (x0, y0) si trovasse nell’origine. A partire dax0 ey0 abbiamo applicatoil metodo delle potenze, ottenendo ogni volta risultati deltipo di quelli riportatinelle figure 1.1 e 1.2, in cui i punti sono 9.

Nella prima figura si vede il poligono iniziale in cui i lati siintrecciano ca-sualmente, poi la matassa si sbroglia e piano piano assistiamo alla nascita dei due“poligoni di accumulazione”P e ¯P entro una trentina di iterazioni: effettivamentela successione da un certo punto in poi oscilla fra questi duepoligoni, i cui verticistanno su un’ellisse.

Ma.......se lasciamo andare ancora avanti il metodo le cosecambiano: come


−0.5 0 0.5−0.5

0

0.5poligono 504

−0.5 0 0.5−1

−0.5

0

0.5

1poligono 558

−0.5 0 0.5−0.5

0

0.5poligono 599

−0.5 0 0.5−0.5

0

0.5poligono 690

Figura 1.2 Esempio 1.5.2: la precisione finita fa convergere i poligoni a un punto

mostra la figura 1.2, i vertici cominciano a spostarsi e pianopiano collassano su unpunto. Cosa sta succedendo? La risposta e semplice: per effetto della precisionefinita la componente dei vettori lungoz1 rispunta fuori e pertanto le due succes-sioni {xk} e {yk} convergono a vettori della formaαz1 eβz1 rispettivamente; laforma diz1 fa il resto.

�

1.5.2 Metodo di Iterazione inversa

Il metodo delle potenze cosı come descritto nel paragrafo precedente non vienemai usato. Viene invece molto utilizzato il metodo di iterazione inversa, notoanche come metodo dellepotenze inverse, che e una variante del metodo dellepotenze originariamente introdotta da H. Wielandt nel 1944nell’articolo [43].

Il metodo si basa sulla seguente osservazione. Siaµ 6= λ1, . . . , λn. Allora lamatrice

B = A− µI

ha autovaloriλ1 − µ, . . . , λn − µ tutti diversi da zero, e gli stessi autovettori diA, mentreB−1 ha autovaloriξi = (λi − µ)−1 peri = 1, . . . , n e ancora gli stessiautovettori diA. L’idea e allora la seguente: seµ e una buona approssimazionedi un autovalore diA, ad esempioλj, e λj e un autovalore semplice, alloraξje l’autovalore dominante diB−1; pertanto il metodo delle potenze applicato aB−1 converge a questo autovalore e contemporaneamente permette di calcolare

30 Capitolo 1

un autovettore ad esso corrispondente7. Piuµ e vicino aλj, piu rapida sara laconvergenza. Una volta che il metodo ci ha portato aξj possiamo ottenere unastima diλj migliore diµ attraverso la formula

λj = µ+1

ξj.

Il metodo di iterazione inversa consiste proprio nell’applicazione del metodo dellepotenze alla matriceB−1. Lo scalareµ in questo contesto e chiamatoshift, perchesi lavora con una matrice il cui spettro e traslato (shifted, in inglese) diµ rispettoallo spettro diA.

La realizzazione del metodo di iterazione inversa richiedeil calcolo del vet-torew = B−1v ad ogni iterazione. Dal punto di vista numerico, questo calcoloviene eseguito risolvendo il sistema lineareBw = v. Sappiamo che in generale larisoluzione di un sistema lineare di dimensionen ha un costo dell’ordine din3. Inrealta, in questo caso i costi possono essere ridotti perche la matriceB = A−µInon cambia nel corso delle iterazioni e pertanto se ne puo calcolare una fattoriz-zazione una volta per tutte prima dell’inizio del ciclo iterativo; ad ogni iterazionesi useranno i fattori per calcolarew. Questo implica un costoO(n3) una sola vol-ta e poiO(n2) ad ogni iterazione. Inoltre, e importante osservare che ilmetododelle potenze inverse viene di solito usato per calcolare autovettori corrispondentiad autovalori di cui si e calcolata una buona stimaµ con altri metodi (ad esempioil metodo QR che vedremo successivamente). In questo contesto, la matriceAviene sempre ridotta inizialmente in forma di Hessenberg, otridiagonale nel casosimmetrico (cfr. paragrafo 1.4). Per matrici di questa forma si possono adattaregli usuali metodi di fattorizzazione, sfruttando la struttura in modo da ridurre icosti; ad esempio in [21] e descritto un algoritmo di Gauss con pivoting parzialeper matrici in forma di Hessenberg superiore che ha un costo dell’ordine di n2,come pure, nello stesso testo, e presentato un algoritmo per la fattorizzazione QRdi una matrice tridiagonale simmetrica con un costo dell’ordine din.

Un’altra osservazione molto interessante dal punto di vista numerico e laseguente. Quanto piuµ e una buona approssimazione di un autovalore diA,tanto piu la matriceB = A − µI e vicina alla singolarita e quindi mal con-dizionata. Molti ricercatori negli anni ’60–’70 ritenevano che questo costituisseun grosso ostacolo all’uso del metodo (nonostante l’evidenza che mostrava buonirisultati) perche e risaputo che l’errore nella risoluzione numerica di sistemi malcondizionati e grande. In realta nel 1979 fu chiarito il mistero in un articolo di G.Peters e J.H. Wilkinson [33], nel quale si dimostra che gli errori nella soluzionetendono a disporsi nella direzione dell’autovettore cercato e scompaiono in virtudella normalizzazione che viene fatta ad ogni iterazione.

In pratica il metodo di Iterazione inversa non viene usato per calcolare gliautovalori, ma e ad oggi il metodo piu utilizzato per calcolare gli autovettori,

7Usandoµ = 0 il metodo convergera all’autovalore diA piu piccolo in modulo, purche questosia isolato rispetto agli altri.


una volta che si siano ottenute buone stime degli autovaloricon qualche altrometodo numerico. Ovviamente, gli algoritmi pratici che realizzano il metodo sonomolto complicati per tenere conto di situazioni non previste nella nostra sempliceintroduzione, ad esempio per riuscire a trattare matrici incui siano presenti gruppidi autovalori molto vicini fra loro e/o matrici difettose, in cui manca un sistemacompleto di autovettori indipendenti. Esempi di tali algoritmi sono le routinesxHSEIN e xSTEIN8, per il caso non simmetrico e simmetrico rispettivamente,della libreria LAPACK. Questi algoritmi sono tali per cui inmolte situazioni esufficiente una iterazione per trovare autovettori molto accurati.

1.5.3 Metodo QR

Poniamoci la domanda seguente. Data una matriceA con autovalori tali che

|λ1| ≥ . . . ≥ |λr| > |λr+1| ≥ . . . ≥ |λn|e possibile generalizzare il metodo delle potenze in modo da calcolare simul-taneamente i primir autovalori? Una generalizzazione di questo tipo dovrebbeprevedere di scegliere una matrice inizialeQ0 ∈ R

n×r ortogonale (l’analogo ma-triciale di un vettore di norma euclidea uguale a 1) che puo essere ottenuta, adesempio, dalla fattorizzazioneQR di una qualsiasi matriceZ ∈ R

n×r, e costruireuna successione di matriciQk ∈ R

n×r ortogonali attraverso due operazioni:

1) premoltiplicazione perA della matrice corrente,2) ortogonalizzazione.

L’idea e schematicamente descritta nel seguente algoritmo, dove l’operazioneal passo 3. del ciclo iterativo generalizza il calcolo dei quozienti di Rayleigh perla stima degli autovalori del metodo delle potenze.

Questo metodo viene chiamato metodo diIterazione simultanea, per met-tere in evidenza il fatto che si itera simultaneamente su pi`u di un vettore, oppureIterazione ortogonaleper evidenziare il fatto che si opera con matrici ortogonali.

Algoritmo 1.5.2. Schema del metodo di Iterazione simultanea1. DataQ0 ∈ R

n×r

2. Perk=1,2, . . .1. PoniZ = AQk−1

2. Calcola la fattorizzazioneZ = QkRk, conQk ∈ Rn×r eRk ∈ R

r×r

3. CalcolaTk = QTkAQk

Perr = 1 il metodo non e altro che il metodo delle potenze, conQk che giocail ruolo di vk eTk ∈ R

r×r quello diσk. Possiamo chiederci a cosa porta questo

8I nomi dei sottoprogrammi di LAPACK sono di solito indicati con l’iniziale “x”, che cor-risponde a “S” per la versione in precisione semplice, “D” per quella in doppia precisione, “C” perquella complessa precisione semplice e infine “Z” per quellacomplessa doppia precisione.

32 Capitolo 1

metodo perr > 1. A questo scopo occorre ricordare il concetto di sottospazioinvariante introdotto nel paragrafo 1.2.1. Si puo infattidimostrare il seguenteteorema:

Teorema 1.5.1.Sia |λ1| ≥ . . . ≥ |λr| > |λr+1| ≥ . . . ≥ |λn| e siaS il sot-tospazio invariante associato a{λ1, . . . , λr}. Allora, seQ0 non e ortogonale aS, la successione{Qk} generata dal metodo di Iterazione simultaneae tale cheRange(Qk), ovvero il sottospazio generato dalle colonne diQk, converge aS,

con la velocita di∣

∣

∣

λr+1

λr

∣

∣

∣

k

.

Per la dimostrazione del teorema rimandiamo a [21, Teorema 7.3.1] per ilcaso non simmetrico e a [21, Teorema 8.2.2] o [40, Teorema 28.2] per il casosimmetrico. Alla luce di questo teorema, considerando il fatto che una base di unsottospazio invariante non e in generale costituita da autovettori, c’e da chieder-si a cosa puo servire il metodo di Iterazione simultanea. Per rispondere, diamoun’occhiata a cosa succede dal punto di vista degli autovalori, ossia delle matriciTk = QT

kAQk. Se la successione{Range(Qk)} tende aS, possiamo intuire chelo spettro diTk tende a{λ1, . . . , λr}. Potremmo allora calcolare questi autovaloricalcolando quelli diT , dove, abusando del termine come abbiamo fatto parlandodel metodo delle potenze, indichiamo conT la matrice a cui leTk “convergono”.In questo modo pero avremmo solo ottenuto di diminuire la dimensione del pro-blema, il che non sarebbe una grande conquista, se non nel caso chen sia moltogrande e si vogliano calcolare soltanto pochi autovalori, ossiar << n.

In realta, il metodo viene utilizzato quandon non e molto grande e si voglionocalcolare tutti gli autovalori, ossiar = n. In questo caso, se si usaQ0 = I la ma-trice T risulta avere una struttura particolare che puo essere sfruttata. In questecondizioni, peraltro, si vede che conviene formulare il metodo in una forma diver-sa, scoprendo che di fatto coincide con ilmetodo QR, uno dei metodi fondamentalidel calcolo numerico, proposto nel 1961 indipendentementedalla ricercatrice so-vietica Vera Nikolaevna Kublanovskaya e dall’inglese JohnG.F. Francis. Questometodo e stato inserito fra i 10 algoritmi “with the greatest influence on the de-velopment and practice of science and engineering in the 20th century”, o anche“one of the jewels on the crown of matrix computations”. Una storia del metodopuo essere trovata nell’articolo [23], e una sua stupenda presentazione in chiavemoderna in [42]. Il metodo QR puo essere schematizzato nel modo seguente:

Algoritmo 1.5.3. Schema del metodo QR1. PoniT0 = A2. Perk=1,2, . . .

1. Calcola la fattorizzazioneTk−1 = QkRk

2. CalcolaTk = RkQk

Prima di tutto stabiliamo l’equivalenza fra questo metodo eil metodo di Itera-zione simultanea. A questo scopo consideriamo innanzitutto l’algoritmo 1.5.2 conr = n. Indichiamo conQk le matrici ortogonali calcolate da questo algoritmo con


Q0 = I, riservando il nomeQk a quelle calcolate dall’algoritmo 1.5.3. Poniamoinoltre T0 = QT

0AQ0 = A e Rk = RkRk−1 . . . R1R0, conR0 = I. Facciamovedere che per ognik le matrici Rk e Qk definiscono una fattorizzazione QR diAk, ovvero

Ak = QkRk. (1.31)

La dimostrazione avviene per induzione. La (1.31) e sicuramente vera perk = 0,essendoA0 = Q0R0 = I. Supponiamo che sia vera perk − 1 e facciamo vedereche vale perk. Infatti per ipotesi induttiva abbiamo

Ak = AAk−1 = AQk−1Rk−1,

e per definizione dell’algoritmo si ha

AQk−1 = QkRk.

Mettendo insieme le due cose, si ottiene

Ak = QkRkRk−1 = QkRk,

ossia la (1.31).Passiamo ora a considerare l’algoritmo 1.5.3 e poniamoQk = Q0Q1 . . . Qk

e Rk = RkRk−1 . . . R1R0 conR0 = Q0 = I. Facciamo vedere, procedendo dinuovo per induzione, che anche in questo caso vale la (1.31) e, inoltre,

Tk = QTkAQk, (1.32)

relazione che vale per definizione nell’algoritmo 1.5.2. Lerelazioni sono en-trambe banalmente vere perk = 0. Supponiamole vere perk−1 e dimostriamoleperk. In questo caso abbiamo, per le due ipotesi induttive,

Ak = AAk−1 = AQk−1Rk−1 = Qk−1Tk−1Rk−1,

e, per definizione nell’algoritmo 1.5.3,

Tk−1 = QkRk,

da cui concludiamo

Ak = Qk−1QkRkRk−1 = QkRk.

Con questo abbiamo dimostrato che vale la (1.31). Per quantoriguarda la (1.32),abbiamo

Tk = RkQk = QTk Tk−1Qk = QT

k QTk−1AQk−1Qk = QT

kAQk,

ovvero la relazione cercata.

34 Capitolo 1

In conclusione, entrambi i metodi 1.5.2 e 1.5.3 producono successioni{Qk},{Rk} e {Tk}, che soddisfano per ognik le due relazioni (1.31) e (1.32). A causadella non unicita della fattorizzazione QR, non possiamo affermare in generaleche le successioni sono le stesse per i due metodi, ma possiamo comunque affer-mare che i due algoritmi sono sostanzialmente la stessa cosa. Il fatto che valga la(1.32) e fondamentale: questa relazione dice infatti che tutte le matriciTk sonoortogonalmente simili adA e pertanto, se la successione ha un limiteT , questamatrice limite ha gli stessi autovalori diA.

L’algoritmo 1.5.3 e la base del metodo QR vero e proprio, chesi ottieneintroducendo tutta una serie di accorgimenti necessari a renderlo se possibile con-vergente (la successione{Tk} deve convergere a unaT con una struttura da cui sipossano evincere facilmente gli autovalori diA) e poi anche efficiente e competi-tivo con gli altri metodi. Fondamentalmente, si deve operare su matrici strutturate(di Hessenberg o tridiagonali) per ridurre il costo delle singole iterazioni, intro-durre degli shifts per forzare e/o accelerare la convergenza e una tecnica di “sgon-fiamento” (o “deflazione”, dall’inglese “deflation”) per ridurre la dimensione delproblema via via che gli autovalori vengono calcolati.

Consideriamo prima di tutto il problema del costo delle iterazioni. Per ridurredrasticamente questo costo e sufficiente trasformare preventivamenteA in formadi Hessenberg, che possiamo supporre, senza perdere in generalita, irriducibile.La fattorizzazione QR di una matrice di Hessenberg, in cui sono da azzeraresoltanto glin − 1 elementi della codiagonale inferiore, puo essere ottenuta inmodo molto efficiente usando trasformazioni di Givens invece che trasformazionidi Householder, che vengono di solito usate per matrici piene e portano ad un cos-toO(n3). Piu precisamente, seT0 e in forma di Hessenberg, la sua fattorizzazioneQR puo essere espressa da

GTn−1 · · ·GT

1 T0 = R1,

doveGj = G(j, j + 1, θj) e la matrice di Givens corrispondente alla scelta(p, q) = (j, j+1) conθj scelto in modo cheGT

j azzeri l’elemento di indici(j, j+1) della matrice a cui e applicata (cfr. algoritmo 1.2.1) e la premoltiplicazione perGT

j puo essere effettuata tramite l’algoritmo 1.2.2.

La matriceQ1 della fattorizzazione QR e data allora daQ1 = G1 · · ·Gn−1.Grazie alla struttura delle matrici di Givens coinvolte, eabbastanza facile vederecheQ1 ha la struttura di Hessenberg superiore e quindi ancheT1 = R1Q1 risultaavere la stessa struttura. In conclusione, la forma di Hessenberg della matriceinizialeT0 viene conservata attraverso le iterazioni.

L’algoritmo 1.5.4 rappresenta uno schema del metodo QR per matrici diHessenberg.


Algoritmo 1.5.4. Schema del metodo QR per matrici di Hessenberg1. TrasformaA in forma di Hessenberg superioreH e poniT = H2. Perk = 1, 2, . . .

1. Perj = 1, . . . , n − 1Calcolacj esj usando l’algoritmo 1.2.1 con datia = tjj e b = tj+1,j

Sovrascrivi aT la matriceGTj T usando l’algoritmo 1.2.2

Fine ciclo suj2. Perj = 1, . . . , n − 1

Sovrascrivi aT la matriceTGj usando l’algoritmo 1.2.3Fine ciclo suj

Fine ciclo suk

Nell’algoritmo abbiamo eliminato gli indici di iterazionee definito una solamatriceT che ad ogni iterazione viene trasformata, passando daTk−1 a Tk. Perogni k, nel primo ciclo suj si definiscono le matriciGj per j = 1, . . . , n − 1e le si applicano via via aT in modo da ottenereRk; nel secondo ciclo suj siapplicano le matriciGj a destra allaT (ovvero allaRk) e si ottiene la nuovaT(ossiaTk). Il costo complessivo di ogni iterazione e dell’ordine din2.

Esempio 1.5.3.Consideriamo la matrice

A =

0.8330 12.674 6.0861 −7.3750−0.2382 6.5607 4.5226 −3.2170−0.2668 2.9621 4.1721 −1.6274−2.3410 9.0454 4.0268 −4.0859

che ha autovalori reali e distinti5.0000, 1.2502, 0.9997 e 0.2299. L’algoritmo1.5.4 applicato allaA converge alla matrice triangolare

T =

5.0000 −14.231 −9.6384 −11.8580 1.2502 −0.3630 −2.40010 0 0.9997 −0.87550 0 0 0.2299

che ha sulla diagonale proprio gli autovalori diA. Constatiamo quindi che l’algo-ritmo ha prodotto la forma di Schur della matriceA, nel senso cheA = QTQT

per una matrice ortogonale Q che potremmo ricostruire usando tutte leQk. �

Esempio 1.5.4.Consideriamo ora una seconda matrice

A =

1 1 1 1 1 11 2 3 4 5 −61 3 6 10 15 211 4 10 20 35 561 5 15 35 70 1261 6 21 56 126 252

36 Capitolo 1

che ha quattro autovalori reali e distinti−1.8341 × 10−2, 0.6260, 16.753, 332.58e due autovalori complessi coniugati0.5301±0.8008i. In questo caso l’algoritmoconverge alla matrice

T =

332.5788 5.2814 −7.1007 5.3581 −0.8875 −0.34400 16.7533 −1.7940 4.9045 −0.2201 0.89110 0 0.9307 1.7914 −0.5132 0.74050 0 −0.4475 0.1295 1.1284 −0.91240 0 0 0 0.6260 −0.16410 0 0 0 0 −0.0183

.

Questa voltaT non e triangolare, e con un semplice ragionamento si capisce cheera impossibile che lo fosse. Ricordiamo infatti che le matrici Tk sono tutte similiallaA e pertanto la matrice limiteT , se triangolare, dovrebbe avere sulla diagonaletutti gli autovalori diA, compresi quelli complessi coniugati. Poiche lavoriamoin aritmetica reale, e impossibile che questo accada. Guardiamo comunque piuda vicino la matriceT . E una matrice triangolare a blocchi con quattro blocchidiagonali di dimensione 1 coincidenti con i quattro autovalori reali diA; il bloccodiagonale2 × 2

(

0.9307 1.7914−0.4475 0.1295

)

ha autovalori complessi coniugati coincidenti con quelli di A. Possiamo pertantodire che l’algoritmo ci ha portato alla forma reale di Schur della matriceA. �

Esempio 1.5.5.Consideriamo ora la matrice

A =

105.7373 −11.2109 40.5759 −16.6085 17.3803 −64.5983363.9555 −41.2617 130.0472 −72.9223 61.3789 −241.1360

−116.5544 16.0762 −35.2613 24.3302 −18.7331 80.3532227.7955 −20.9683 85.3947 −46.8571 41.1281 −142.4106

−106.9441 11.8851 −38.0521 21.1965 −16.7694 71.3905−115.3461 8.8529 −45.4448 28.0879 −21.0790 80.4122

i cui autovalori sono:46.9995,−11.9997, 10.0002, 1.0000 e±5.2345. L’algorit-mo converge alla matrice triangolare a blocchi

T =

46.9995 441.8497 145.6208 259.9014 −120.0433 288.32720 −11.9997 −1.1263 −5.9210 4.5155 −10.10580 0 10.0002 −3.3827 −4.6832 −2.73020 0 0 3.1515 −8.8030 −2.83830 0 0 −1.9843 −3.1515 −0.15170 0 0 0 0 1.0000

.

I blocchi diagonali1 × 1 coincidono con i primi quattro autovalori diA, mentregli autovalori del blocco2 × 2

(

3.1515 −8.8030−1.9843 −3.1515

)


sono±5.2345. Pertanto, l’algoritmo conduce a una matrice da cui si calcolanoagevolmente tutti gli autovalori diA, ma non alla forma di Schur che in questocaso dovrebbe presentare unaT triangolare superiore. �

Esempio 1.5.6.Consideriamo infine la matrice di Hessenberg superiore

A =

(

0 0 11 0 00 1 0

)

i cui autovalori sono1.0000 e −0.5000 ± 0.8660i. Applicando l’algoritmo siottiene una successione{Tk} non convergente, con

T3i =

(

0 0 11 0 00 1 0

)

, T3i+1 =

(

0 0 −11 0 00 −1 0

)

, T3i+2 =

(

0 0 −1−1 0 0

0 1 0

)

peri = 0, 1, 2, . . .. �

Alla luce degli esempi precedenti, possiamo chiederci quali siano le proprietadi convergenza dell’algoritmo 1.5.4. Il risultato dell’esempio 1.5.3 si spiega conil seguente teorema 1.5.2, valido nel caso in cui gli autovalori siano tutti reali edistinti. Il teorema vale piu in generale per l’algoritmo 1.5.3 e non solo per matricireali, ma anche per matrici a elementi complessi [5]. L’ipotesi sull’esistenza dellafattorizzazione LU diX−1 e un’ipotesi tecnica, che potrebbe essere sostituita daaltre ipotesi a seconda della strada seguita per la dimostrazione. Le matrici di fasedi cui si parla nell’enunciato sono matrici diagonali con elementi diagonali ugualia±1 e intervengono per tenere conto della non unicita della fattorizzazione QR(cfr. [5] o [21]).

Teorema 1.5.2.Supponiamo cheA ammetta autovalori reali e distinti, con

|λ1| > |λ2| > . . . > |λn|, (1.33)

e siaX−1AX = Λ la sua diagonalizzazione. Supponiamo cheX−1 ammettafattorizzazione LU. Allora, esistono delle matrici di faseSk, conk ≥ 0, tali che

limk→∞

SkTkSk−1 = T

conT triangolare superiore con gli autovalori diA sulla diagonale, e

limk→∞

Sk−1QkSk = I.

In sostanza il teorema dice che quando gli autovalori sono tutti distinti inmodulo il metodo QR converge alla forma di Schur diA. Per quanto riguardala velocita di convergenza, si puo dimostrare che il metodo converge linearmente

38 Capitolo 1

con velocitaO(

mini

∣

∣

∣

λi+1

λi

∣

∣

∣

)

e quindi puo essere lento quando ci sono autovalori

vicini in valore assoluto. Questo non ci sorprende, trattandosi in fondo in fondodi un’evoluzione del metodo delle potenze.

Risultati di convergenza nel caso che l’ipotesi (1.33) non sia verificata si pos-sono trovare gia nel libro di Wilkinson [45]. Ad esempio, seesiste un certo nu-mero di autovalori reali coincidenti, allora si puo ancoradimostrare la convergen-za alla forma di Schur. D’altra parte, l’esempio 1.5.5 ci mostra che questo none piu vero in presenza di autovalori reali opposti. Per quanto riguarda il caso incui siano presenti autovalori complessi coniugati, abbiamo due esempi contrad-dittori: nell’esempio 1.5.4 il metodo ci porta alla forma diSchur reale, mentrenell’esempio 1.5.6 il metodo non converge.

Non vale la pena investigare ulteriormente questo punto perche il metodo QRviene realizzato di fatto in una forma diversa da quella semplice vista finora: perforzare e/o accelerare la convergenza si utilizza una strategia dishifting, la stessaidea che abbiamo visto nel metodo di iterazione inversa. Introduciamo questenuove forme del metodo, e poi ci occuperemo delle proprietadi convergenza. Inprima istanza consideriamo il seguente schema delmetodo QR con shift, nel qualeindichiamo conσk lo shift scelto ad ogni iterazione.

Algoritmo 1.5.5. Schema del metodo QR con shift1. TrasformaA in forma di Hessenberg superioreH e poniT0 = H2. Perk=1,2, . . .

1. Scegliσk “vicino ad un autovalore di A”2. Calcola la fattorizzazioneTk−1 − σkI = QkRk

3. CalcolaTk = RkQk + σkIFacciamo vedere prima di tutto che le matriciTk sono ancora tutte ortogonal-

mente simili fra loro e quindi adA. Infatti

Tk = RkQk + σkI = QTk (QkRk + σkI)Qk = QkTk−1Qk. (1.34)

Inoltre, poiche la matriceTk−1 − σkI e ancora di Hessenberg, le considerazionifatte per il metodo senza shift ci dicono che tutte le matriciTk sono in forma diHessenberg superiore.

Il problema e come scegliere lo shiftσk. Nelle considerazioni che seguonoeliminiamo gli indici di iterazione e indichiamo conT e T la vecchia e la nuovamatrice. Seσ e un autovalore diA, alloraT − σI e singolare e altrettanto valeperR, che quindi ha un elemento diagonale uguale a zero. Supponiamo che siarnn = 0. Allora, la rigan-esima diRQ e zero e di conseguenza la rigan-esimadi T e σeTn , doveen e l’n-esimo vettore della base canonica.T e pertanto unamatrice di Hessenberg non irriducibile in cui e stato isolato l’autovaloreσ, e ilmetodo puo proseguire lavorando sul minore principale diT di dimensionen−1.

Seσ non e esattamente un autovalore diA, ma e comunque (molto) vicinoad uno di essi, diciamoσ ≃ λj, allora ci possiamo aspettare che la matriceT siaquasi non irriducibile. Per continuare l’esempio precedente, ci possiamo aspettare


che tn,n−1 ≃ 0 (e tnn ≃ λj). Se al procedere delle iterazioniσ viene sceltosempre piu vicino aλj, ci possiamo aspettare chetn,n−1 tenda a zero. Questaaspettativa puo basarsi sul seguente ragionamento, da cuiemerge che il metodoQR sta facendo un’iterazione del metodo di Iterazione inversa. Infatti, indichiamoconq l’ultima colonna diQ e supponiamo cheσ sia un autovalore diA ernn = 0.Allora si ha

qTT = qT (QR+ σI) = eTnR+ σqT = σqT

e pertantoq e un autovettore sinistro diT . Supponiamo ora di nuovo cheσ non siaesattamente un autovalore diAma vicino ad uno di essi. Allora, daT−σI = QRsegue

(T − σI)R−1 = Q ⇒ R(T − σI)−1 = QT

⇒ (T T − σI)−1RT = Q⇒ (T T − σI)−1rnnen = q,

ovvero q e quello che otterremmo applicando un’iterazione del metodo di Itera-zione inversa aT T con shiftσ, e pertantoq e vicino ad un autovettore sinistro diT , tanto piu quanto piuσ e una buona stima dell’autovalore.

Come scegliereσ perche le aspettative si realizzino? La prima idea e quella diprendere ad ogni iterazioneσ = tnn, in modo da forzare la convergenza diTk aduna matriceT non irriducibile con un autovalore diA in posizione(n, n). Questastrategia funziona bene in generale se la matriceA ammette qualche autovalorereale isolato; si puo anche far vedere che l’elemento in posizione(n, n− 1) di Tk

tende a zero con velocita quadratica (vedi [8, pag. 162] e [21, pag. 354]).

L’algoritmo 1.5.5 puo creare difficolta seA ha una coppia di autovalori com-plessi coniugati. Infatti in questo caso puo accadere che la sottomatrice

D =

(

tn−1,n−1 tn−1,n

tn,n−1 tnn

)

(1.35)

abbia autovalori complessi e quindi prendere come shifttnn (che e reale) puo nonessere una buona strategia. Per questo motivo gia Francis nel suo articolo del1961 propose di usare una strategia piu complessa, dettadoppio shift implicito,nella quale si calcola la matriceTk a partire daTk−1 mediante due passi di shiftconsecutivi, usando come shifts gli autovaloriα1 e α2 di D. Vedremo che que-sta strategia puo essere realizzata in aritmetica reale anche quando siano presentiautovalori complessi. Il motivo dell’aggettivo “implicito” sara chiarito piu avanti.

Indichiamo conT (0) e T (2) le matrici Tk−1 e Tk rispettivamente. I duepassaggi di shift sono espressi da

T (0) − α1I = Q1R1

T (1) = R1Q1 + α1IT (1) − α2I = Q2R2

T (2) = R2Q2 + α2I(1.36)

40 Capitolo 1

Procedendo come in (1.34), si deduce cheT (1) = QT1 T

(0)Q1 eT (2) = QT2 T

(2)Q2,da cui segue

T (2) = QT2Q

T1 T

(0)Q1Q2. (1.37)

Vedremo che l’algoritmo del QR con doppio shift implicito sibasa sul calcolo“furbo” della matriceZ = Q1Q2. Prima di vedere come funziona l’algoritmo,chiediamoci se le operazioni (1.36) ci fanno veramente mantenere l’aritmeticareale anche seα1 eα2 sono complessi coniugati. La risposta e positiva. Infatti,postoR = R2R1, si ottiene

ZR = (Q1Q2)R2R1 = Q1(T(1) − α2I)R1

= Q1(R1Q1 + α1I − α2I)R1

= Q1R1Q1R1 + (α1 − α2)Q1R1

= (T (0) − α1I)(T (0) − α1I) + (α1 − α2)(T(0) − α1I)

= (T (0) − α1I)[T (0) − α1I + (α1 − α2)I]= (T (0) − α1I)(T (0) − α2I)= (T (0))2 − (α1 + α2)T

(0) + α1α2I ≡M.

La matriceM e reale perche i due coefficientiα1 +α2 eα1α2 sono entrambi realianche quandoα1 eα2 sono complessi coniugati. Pertanto e possibile scegliereQ1

eQ2 in modo cheZ edR siano reali; con questa scelta ancheT (2) e reale perche

T (2) = ZTT (0)Z. (1.38)

In pratica potremmo allora mantenere l’aritmetica reale evitando di usare leformule (1.36) e seguendo invece la strada alternativa che consiste nel costruireM , calcolarne la fattorizzazione QR

M = ZR (1.39)

e calcolare infineT (2) con la (1.38). Il calcolo della fattorizzazione (1.39) ha uncosto computazionaleO(n3). Questo costo puo essere ridotto aO(n2) seguendoun approccio leggermente diverso che si basa su un importante teorema di algebralineare, noto comeTeorema del Q implicito. Esistono molte versioni di questoteorema; a noi interessa quella relativa alle matrici di Hessenberg irriducibili.

Teorema 1.5.3.SiaA = QHQT e A = V GV T , conH e G di Hessenbergsuperiore irriducibili eQ = (q1 . . . qn) eV = (v1 . . . vn) ortogonali. Seq1 = v1,allora qi = ±vi per i = 2, . . . , n.

Dimostrazione. PoniamoW = V TQ = (w1 · · ·wn). Dall’ortogonalita diQ e V e dall’ipotesiq1 = v1 segue chew1 = e1; inoltre si vede facilmente cheGW = WH, da cui segue

Gwi−1 =

i∑

j=1

hj,i−1wj ,


ovvero

hi,i−1wi = Gwi−1 −i−1∑

j=1

hj,i−1wj, (1.40)

per i = 2, . . . , n. Peri = 2 la relazione (1.40) dice cheh21w2 = Gw1 − h11w1;essendoG di Hessenberg superiore ew1 = e1, segue chew2 ha soltanto i primidue elementi diversi da zero. Per induzione, si arriva a dimostrare cheW e trian-golare superiore. PoicheW e anche ortogonale, essa deve necessariamente esserediagonale conwi = ±ei peri = 2, . . . , n, il che equivale a dire cheqi = ±vi peri = 2, . . . , n.

Avendo in mente questo teorema, si puo costruire a costoO(n2) una matriceZ che soddisfa contemporaneamente (1.39) e (1.38); i passi daeseguire sono iseguenti:

• Si calcola la prima colonna diM ; indicando contij gli elementi diT (0), essa edata da

Me1 =

t211 + t12t21 − (α1 + α2)t11 + α1α2

t21(t11 + t22 − (α1 + α2))t21t32

0...0

• Come se si volesse calcolare la fattorizzazione (1.39), si individua una matricedi HouseholderP1 tale cheP1Me1 sia un multiplo die1. Il costo dell’ operazioneeO(1) perche si tratta di annullare soltanto due elementi di un vettore inR

n.

•Si calcolano delle matrici di HouseholderP2, . . . , Pn−1 tali che la prima colonnadi Z1 = P1 · · ·Pn−1 resti uguale alla prima colonna diP1 (ovvero alla primacolonna diZ) e la matriceZT

1 T(0)Z1 sia di Hessenberg superiore. A questo punto

il teorema 1.5.3 implica cheZ1 coincide conZ a meno dei segni delle colonnesuccessive alla prima e quindiZT

1 T(0)Z1 coincide conT (2).

Come costruireZ1 rispettando l’impegno a mantenere il costo computazionaledell’ordine din2? L’algoritmo che viene utilizzato a questo scopo e un algoritmocosiddettobulge-chasing(che lima le sporgenze), gia proposto da Francis nell’ar-ticolo del 1961, che funziona nel modo che andiamo a descrivere. Osserviamoprima di tutto che la matriceP1 ha la forma

P1 =

(

P1 00 In−3

)

,

42 Capitolo 1

con P1 matrice di Householder, e pertanto calcolandoP T1 T

(0)P1 si modificanosoltanto le prime 3 righe e colonne diT (0); in pratica, vengono introdotti deglielementi diversi da zero nelle posizioni (3,1), (4,1) e (4,2), come evidenziato nellafigura sottostante (fatta per il cason = 6):

P T1 T

(0)P1 =

∗ ∗ ∗ ∗ ∗ ∗∗ ∗ ∗ ∗ ∗ ∗⋄ ∗ ∗ ∗ ∗ ∗⋄ ⋄ ∗ ∗ ∗ ∗0 0 0 ∗ ∗ ∗0 0 0 0 ∗ ∗

.

Questa prima trasformazione ha di fatto introdotto una sporgenza nella forma diHessenberg. Le trasformazioni successive sposteranno la sporgenza via via piuin basso fino a farla scomparire. A questo scopo scegliamoP2 come una trasfor-mazione di Householder che annulli gli elementi in posizione (3,1) e (4,1); questaP2 coincide con la matrice identita al di fuori delle righe e colonne 2, 3 e 4.Avremo quindi

P T2 P

T1 T

(0)P1P2 =

∗ ∗ ∗ ∗ ∗ ∗∗ ∗ ∗ ∗ ∗ ∗0 ∗ ∗ ∗ ∗ ∗0 ⋄ ∗ ∗ ∗ ∗0 ⋄ ⋄ ∗ ∗ ∗0 0 0 0 ∗ ∗

.

A questo punto,P3 deve azzerare le posizioni (4,2) e (5,2) ed e quindi l’identitaal di fuori delle righe 3, 4 e 5, ossia

P T3 P

T2 P

T1 T

(0)P1P2P3 =

∗ ∗ ∗ ∗ ∗ ∗∗ ∗ ∗ ∗ ∗ ∗0 ∗ ∗ ∗ ∗ ∗0 0 ∗ ∗ ∗ ∗0 0 ⋄ ∗ ∗ ∗0 0 ⋄ ⋄ ∗ ∗

.

Continuando nello stesso modo, scegliamoP4 in modo da annullare le posizioni(5,3) e (6,3).P4 coincide con l’identita al di fuori delle righe 4, 5 e 6, da cui segue

P T4 P

T3 P

T2 P

T1 T

(0)P1P2P3P4 =

∗ ∗ ∗ ∗ ∗ ∗∗ ∗ ∗ ∗ ∗ ∗0 ∗ ∗ ∗ ∗ ∗0 0 ∗ ∗ ∗ ∗0 0 0 ∗ ∗ ∗0 0 0 ⋄ ∗ ∗

.

La figura mostra che nel cason = 6 siamo giunti dopo 4 trasformazioni a doverannullare soltanto un elemento nell’ultima riga. Pern generico si arriva a questa


situazione dopon − 2 trasformazioni. L’ultima trasformazione puo essere sceltacome matrice di Givens

Pn−1 =

( In−2

c s−s c

)

conc eds tali da annullare l’elemento in posizione(n, n− 2). In questo modo lamatrice

T (2) = P Tn−1 · · ·P T

1 T(0)P1 · · ·Pn−1

risulta di Hessenberg superiore.Osserviamo adesso che la prima colonna diZ1 = P1 · · ·Pn−1 coincide con

la prima colonna diP1 perche le trasformazioni successive non la modificano.Questo potremmo dirlo anche della matriceZ costruita nel processo di fattoriz-zazione QR diM a partire daP1. Abbiamo cosı ottenuto lo scopo desiderato.

Possiamo notare che l’algoritmo ora tratteggiato permettedi fare l’iterazionedi QR con doppio shift senza bisogno di formare esplicitamente le due matriciT (0) − α1I e T (1) − α2I; questo e il motivo per cui si usa l’aggettivo implicitonella denominazione dell’algoritmo.

Il metodo QR nella forma appena descritta converge in generale alla forma diSchur reale diAmolto velocemente, dimostrando con cio che la scelta degliauto-valori della sottomatrice (1.35) come shifts funziona fin dall’inizio, anche quandoniente assicura che essi siano vicini ad autovalori diA. L’esperienza dimostrache tipicamente l’intero spettro diA viene calcolato con un numero di iterazionidell’ordine din. Siccome ogni iterazione costaO(n2), qualcuno vede il metodocome un metodo “diretto” con un costo complessivoO(n3).

D’altra parte, questo e il classico esempio di metodo numerico universalmenteutilizzato perche le sue prestazioni sono ottime, manon esistono risultati teori-ci che ne garantiscano la convergenza per qualsiasi matrice. Al contrario, sononote matrici su cui il metodo non funziona, come quella dell’esempio 1.5.6 chenon risulta modificata dal metodo percheα1 eα2 sono entrambi uguali a zero. Persopperire alle situazioni in cui il metodo non converge o cicla, sono state suggeritestrategie di doppio shift alternative che sono inglobate nel migliore software perirrobustire l’algoritmo. Ad esempio, una strategia usata `e il cosiddettoshift diWilkinson[29]: se dopo 10 iterazioni con il doppio shift di Francis ancora nes-sun elemento sottodiagonale e diventato talmente piccoloda essere trascurabile,ovvero il metodo non si e attestato su un autovalore (o una coppia di autovalori) diA, gli shifts usualiα1 eα2 sono sostituiti per una volta da due valori che siano inqualche modo legati alla norma della matrice; in particolare, Wilkinson suggeriscedi usareα1 eα2 tali che

α1 + α2 = 1.5(|tn,n−1| + |tn−1,n−2|) , α1α2 = (|tn,n−1| + |tn−1,n−2|)2.

Dopo altre 10 iterazioni, ci si ripone il problema e si fa eventualmente un altro shifteccezionale. Se dopo 30 iterazioni ancora non si e ottenutolo scopo, ci si arrende

44 Capitolo 1

e si dichiara non convergenza. Lo shift di Wilkinson permette di recuperare laconvergenza in molti casi, ma non in tutti. Ad esempio, consideriamo matricidella forma

0 1 0 01 0 −δ 00 δ 0 10 0 1 0

con δ piccolo, ma non troppo vicino alla precisione di macchina (ad esempio,δ = 10−8 lavorando in doppia precisione). Gli autovalori esatti sono complessiconiugati a coppie

±√

1 − δ2

4± i

δ

2.

Essendoδ piccolo, questi autovalori sono vicini a una coppia di autovalori re-ali doppi. Come per l’esempio 1.5.6, il doppio shift di Francis non funziona.Ma purtroppo lavorando in aritmetica finita non funziona neanche lo shift diWilkinson, e si sono studiate altre possibili scelte ad hoc.

Molte altre cose sarebbero da dire riguardo alla realizzazione effettiva diquesto metodo, in particolare riguardo ai criteri adottatiper decidere che unamatrice di Hessenberg ha qualche elemento codiagonale “uguale a zero” e chequindi il problema puo essere “sgonfiato”; per questioni dispazio e di semplicitaci limitiamo a rimandare i lettori interessati a [21] e [42].

1.6 Metodi per matrici simmetriche

Per calcolare gli autovalori di una matrice simmetrica, si possono usare ovvia-mente i metodi per matrici generali. Abbiamo gia visto, ad esempio, che il meto-do di Iterazione inversa puo essere molto efficiente nel caso di matrici tridiagonalisimmetriche. Nel prossimo paragrafo vedremo che anche il metodo QR puo es-sere realizzato in modo tale da trarre vantaggio della simmetria della matrice, inparticolare se la matrice viene preventivamente trasformata in forma tridiagonale.

D’altra parte, esistono ad oggi metodi molto efficienti pensati per il calcolodegli autovalori di matrici di questo tipo. Un panorama aggiornato e chiarissimodi questi metodi si puo trovare nell’articolo [11], nel quale si fa il punto dellasituazione tramite un confronto fra i risolutori presenti nella libreria LAPACK.

I metodi presi in considerazione nell’articolo sono quattro: il metodo QRper problemi simmetrici, il metodo di Bisezione, il metodo Divide–and–Conquer,e il metodo MRRR (Multiple Relatively Robust Representations). Ciascuno diquesti puo essere preferibile agli altri per alcune tipologie di matrici o in alcuniambienti di calcolo. Noi descriveremo i metodi di Bisezionee Divide and Con-quer e tralasceremo il metodo MRRR che richiede basi di algebra lineare e anal-isi numerica che vanno oltre a quelle previste per il potenziale uditorio di questiappunti.

Oltre a questi metodi introdurremo anche due metodi “antichi”. Il primo enoto come metodo dellaIterazione quoziente di Rayleighe mette insieme le idee


del metodo delle potenze e dell’Iterazione inversa ottenendo un metodo a conver-genza cubica. Anche se oggigiorno non e piu utilizzato, resta un “bel pezzo dimatematica”. Il secondo e ilmetodo di Jacobi, l’unico metodo che da come risul-tato tutti gli autovalori e tutti gli autovettori; esso non puo competere in generalecon i metodi piu moderni dal punto di vista dell’efficienza,ma in talune situazioniviene ancor oggi utilizzato.

1.6.1 Metodo QR per matrici simmetriche tridiagonali

SeA e simmetrica, non e necessario considerare strategie di doppio shift per ilmetodo QR, non essendoci il problema dell’esistenza di eventuali autovalori com-plessi coniugati. Consideriamo allora l’algoritmo di base1.5.5 del metodo QR conshift e osserviamo prima di tutto che seT0 e tridiagonale simmetrica, anche tuttele matriciTk generate dall’algoritmo mantengono questa struttura. La scelta diσad ogni iterazione potrebbe essere quella gia suggerita nel paragrafo precedente,σ = tnn. In realta, la strategia suggerita da J.H. Wilkinson [46] euniversalmenteaccettata e un po’ diversa. SeT ha la forma

T =

a1 b1 0 0 · · · 0b1 a2 b2 0 · · · 00 b2 a3 b3 · · · 0

.. . . . . .. .. . . .. . bn−1

0 0 0 0 bn−1 an

,

si prendeσ uguale all’autovalore della sottomatrice

D =

(

an−1 bn−1

bn−1 an

)

piu vicino aan; con calcoli abbastanza semplici si vede che, postod = (an−1+an)2 ,

questo autovalore e dato da

σ = an + d− sign(d)√

d2 + b2n−1.

Come abbiamo visto nel caso delle matrici non simmetriche, anche in questocaso e possibile realizzare un’iterazione del metodo senza formare esplicitamentela matriceT − σI. Infatti si riesce, basandosi sul teorema del Q implicito, acostruire un algoritmo di tipobulge-chasing, in cui il costo dell’aggiornamentodella matriceT diventaO(n), come mostrato in [21]. Vediamo per sommi capicome si passa dalla matrice attualeT alla successiva matriceT = RQ + σI.Data la struttura diT conviene utilizzare trasformazioni di Givens invece che diHouseholder.

46 Capitolo 1

Sianoc = cosθ1 es = sinθ1 tali che

(

c s−s c

)T (a1 − σb1

)

=

(

∗0

)

e poniamoG1 = G(1, 2, θ1) ∈ Rn×n, allora e facile vedere che la matrice

GT1 TG1 ha la forma “quasi-tridiagonale” evidenziata nella figura sottostante (per

il cason = 6)

GT1 TG1 =

∗ ∗ ⋄ 0 0 0∗ ∗ ∗ 0 0 0⋄ ∗ ∗ ∗ 0 00 0 ∗ ∗ ∗ 00 0 0 ∗ ∗ ∗0 0 0 0 ∗ ∗

,

dove gli elementi diversi da zero in posizione (1,3) e (3,1) nascono a causa dellapremoltiplicazione perGT

1 e della postmoltiplicazione perG1 rispettivamente.E inoltre immediato verificare che la prima colonna diG1 e uguale alla primacolonna della matriceQ che otterremmo facendo la fattorizzazione QR diT −σI.

Per poter utilizzare il teorema del Q implicito, dobbiamo essere in grado ditrovare altren − 2 matrici di GivensG2, . . . , Gn−1 tali cheZ = G1G2 · · ·Gn−1

abbia la prima colonna uguale a quella diG1 (e quindi a quella diQ) eZTTZ siatridiagonale. Il primo scopo e sicuramente ottenuto seGi = G(i, i + 1, θi) peri = 2, . . . , n − 1; inoltre, scegliendo in modo opportuno i valoriθi possiamo farscorrere il “⋄” verso il basso in modo da farlo sparire daZTTZ, ottenendo cosıil secondo scopo. A questo punto, in virtu del teorema del Q implicito, ZTTZ eproprio laT desiderata.

Le proprieta di convergenza del metodo QR descritte nel paragrafo prece-dente per matrici non simmetriche continuano a valere anchenel caso simmetrico;in questo caso, ovviamente, se c’e convergenza, la matricelimite e diagonale. Dinuovo, il buon funzionamento del metodo dipende fortementeda alcuni dettaglipratici per i quali rimandiamo a [21].

1.6.2 Metodo Iterazione quoziente di Rayleigh

Il metodo di cui parliamo in questo paragrafo e intimamentelegato al metododi Iterazione inversa e al metodo QR. Ricordiamo che nel metodo di Iterazioneinversa si applica il metodo delle potenze alla matrice(A−µI)−1 ai fini di calco-lare l’autovalore diA piu vicino aµ e, soprattutto, l’autovettore ad esso associato.Lo shift µ e scelto a priori e mantenuto fisso per tutte le iterazioni. L’idea delmetodo Iterazione quoziente di Rayleigh e quella di cambiareµ ad ogni iterazio-ne (come avviene nel metodo QR), usando ogni volta lo shift uguale all’ultimaapprossimazione dell’autovalore calcolata. Il metodo e dovuto essenzialmente aA.M. Ostrowski [31], il quale ha scritto diversi articoli sull’argomento negli anni1957–1959.


Schematicamente, un’iterazione del metodo puo essere descritta dalle for-mule seguenti:

1. w = (A− σk−1I)−1vk−1

2. vk = w‖w‖

3. σk = vTk Avk

(1.41)

Evidentemente, al primo passo occorre risolvere un sistemalineare; poichela matrice cambia conk, la fattorizzazione non puo essere fatta una volta per tutteprima dell’inizio delle iterazioni, come invece si puo fare nel metodo di Iterazioneinversa. D’altra parte, la velocita di convergenza di questo metodo e molto elevata.Vale infatti il seguente risultato

Teorema 1.6.1.SiaA una matrice simmetrica. Allora:

i) Il metodo Iterazione quoziente di Rayleigh converge a unacoppia autovalore–autovettore quasi per ogniv0;

ii) Se il metodo converge, la velocita e cubica, ovvero

‖vk+1 − (±x)‖ = O(‖vk − (±x)‖3)

e|σk+1 − (λ)‖ = O(|σk − λ)|3),

doveλ ex sono l’autovalore e l’autovettore a cui si converge.

La dimostrazione del teorema e molto complicata. Rimandiamo chi fosseinteressato al libro [32] e ci limitiamo, come in [21], a far vedere che le tesi sonovere in un caso molto particolare.

SiaA =

(

λ1 00 λ2

)

, con λ1 > λ2. Gli autovettori sonox1 = (1, 0)T

e x2 = (0, 1)T . Sia vk ∈ R2 un vettore di norma euclidea uguale a 1; allora

possiamo scrivere

vk =

(

cksk

)

,

con c2k + s2k = 1. Partendo davk, calcoliamo primaσk, poi w e infinevk+1.Facendo i calcoli, si ottieneσk = λ1c

2k + λ2s

2k e

w =

ck

λ1−σk

sk

λ2−σk

=

ck

λ1(1−c2k)−λ2s2

k

sk

λ2(1−s2k)−λ1c2

k

=1

λ1 − λ2

ck

s2k

−sk

c2k

,

con

‖w‖ =1

λ1 − λ2

√

c6k + s6k

s2kc2k

.

48 Capitolo 1

Normalizzandow si ottiene allora

vk+1 =

(

ck+1

sk+1

)

=

c3k√

c6k+s6

k

−s3k√

c6k+s6

k

,

da cui deduciamo che la successione converge cubicamente ax1 oppurex2 purchesia|c0| 6= |s0|.

Il problema di questo metodo e che converge ad un autovalore, ma non si saa priori a quale. Allora in linea di principio puo essere usato per calcolare tuttigli autovalori, purche accompagnato da un processo di sgonfiamento del proble-ma. Un modo puo essere il seguente. Dato un autovaloreλ e un corrispondenteautovettorex con‖x‖2 = 1, si costruisce una matrice ortogonaleP la cui primacolonna coincide conx 9. SiaP = [x P1]. Allora

P TAP =

(

xTAx xTP1

P T1 Ax P1TAP1

)

=

(

λ 0λP T

1 x A1

)

=

(

λ 00 A1

)

.

A questo punto si puo lavorare sulla matriceA1. Infatti, gli autovalori diA1 sono

autovalori diA; inoltre, sey ∈ Rn−1 e autovettore diA1, il vettorev = P

(

0y

)

e autovettore diA.Il metodo Iterazione quoziente di Rayleigh e anche strettamente legato al

metodo QR con shift. Si dimostra infatti che se si usa lo shiftσk = tnn nel metodoQR e si inizializza il metodo Iterazione quoziente di Rayleigh con il vettoreen,allora le successioni{σk} generate dai due metodi sono identiche [8].

1.6.3 Metodo Divide–and–conquer

Il nome di questo metodo e la traduzione inglese del famoso motto latino “Divideet impera”. Il metodo si applica a matrici tridiagonali simmetriche irriducibili del-la forma (1.22) e fu introdotto nel 1981 da Cuppen [7]. L’ideae quella di scom-porre il problema in due sottoproblemi di dimensione pari a (circa) meta di quelladel problema originale, calcolare autovalori e autovettori di questi sottoproblemi,e risalire da questi ai risultati per il problema originale.L’idea si puo applicare inmodo ricorsivo, in linea di principio fino a ricondursi a problemi di dimensione 2o 1 per i quali si possono fare i calcoli in modo diretto. In LAPACK si ferma l’ap-plicazione ricorsiva a sottoproblemi di una dimensione sufficientemente piccolada rendere preferibile risolverli con il metodo QR invece che continuare con il

9Ad esempio,P puo essere presa come la matrice di Householder tale chePx = e1. Allora,x = P T e1 e, essendoP simmetrica,Pe1 = x.


divide–and–conquer. Questa dimensione dipende dalla macchina su cui si sta la-vorando, e tiene conto sia dell’occupazione di memoria che del tempo di calcolo.Il procedimento e adattissimo alla realizzazione su elaboratori paralleli, come pro-posto in [13] e riassunto in [21]. Questa caratteristica ha portato ad approfondirelo studio delle proprieta numeriche del metodo (convergenza, accuratezza) e aproporre varianti con ottime proprieta ([3], [34]).

Per semplicita, consideriamon pari, diciamon = 2m, e siaT la matricetridiagonale simmetrica di cui vogliamo calcolare autovalori e autovettori. In-dichiamo cone1 ed em il primo e l’ultimo vettore della base canonica inRm, econv il vettore

v =

(

eme1

)

.

Allora possiamo scrivere

T =

(

T1 βmemeT1

βme1eTm T2

)

=

(

T1 0

0 T2

)

+ βmvvT ,

dove

T1 =

α1 β1 0 · · · 0β1 α2 β2 · · · 0

. . ... .

. . ... . . . . βm−1

0 0 0 βm−1 αm

e

T2 =

αm+1 βm+1 0 · · · 0βm+1 αm+2 βm+2 · · · 0

. . . . .. . . .. .. . . . βn−1

0 0 0 βn−1 αn

,

conαm = αm −βm e αm+1 = αm+1 −βm. Piu in generale, possiamo introdurreun parametroθ con cui definire

v =

(

emθe1

)

e ottenere

T =

(

T1 0

0 T2

)

+βm

θvvT ,

con αm = αm − βm

θe αm+1 = αm+1 − θβm. Il parametroθ puo essere scelto

in modo da evitare potenziali difficolta numeriche legate al fenomeno della can-cellazione nel calcolo diαm e αm+1. Seαm e αm+1 hanno lo stesso segno, siprendeθ = ±1 in modo che nel calcolo diαm e αm+1 non ci siano sottrazioni.

50 Capitolo 1

Altrimenti, si sceglie il segno diθ in modo da evitare sottrazione in uno dei duecalcoli e il suo valore assoluto in modo da evitare sottrazione fra numeri troppovicini nell’altro.

Ora supponiamo di calcolare le decomposizioni di SchurT1 = Q1D1QT1 e

T2 = Q2D2QT2 . Allora, postoρ = βm

θ, possiamo scrivere

T =

(

Q1D1QT1 0

0 Q2D2QT2

)

+ ρvvT

=

(

Q1 00 Q2

)(

D1 00 D2

)(

QT1 0

0 QT2

)

+ ρvvT .

Ponendo

Q =

(

Q1 00 Q2

)

e D =

(

D1 00 D2

)

,

e sviluppando la relazione appena scritta, si ottiene

T = QDQT + ρvvT = QDQT + ρQQT vvTQQT

= QDQT + ρQzzTQT ,

conz = QT v, ovvero

T = Q(D + ρzzT )QT . (1.42)

Vediamo pertanto che, noti autovalori e autovettori diT1 e T2, possiamo otteneregli autovettori diT senza effettuare alcun calcolo e i suoi autovalori calcolandogli autovalori diD + ρzzT . Questa e una matrice speciale in quanto somma diuna matrice diagonale e una matrice di rango 1 (osserviamo che ρ 6= 0 a causadell’irriducibilita di T ) e i suoi autovalori possono essere ottenuti con poco sforzocomputazionale tenendo conto del seguente teorema.

Teorema 1.6.2.SiaV T (D + ρzzT )V = diag(µ1, . . . µn) la decomposizione diSchur diD + ρzzT . Supponiamoz1, z2, . . . , zn 6= 0 eD = diag(d1, . . . , dn) cond1 > d2 > . . . dn. Allora:

i) µ1, . . . µn sono soluzioni dell’equazionef(µ) = 0, dove

f(µ) = (1 + zT (D − µI)−1z)1 + ρn∑

i=1

z2i

di − µ.

ii) Seρ > 0, allora µ1 > d1 > µ2 > d2 > . . . > dn, altrimentid1 > µ1 > d2 >µ2 > . . . > µn;

iii) La colonna i-esima diV e multiplo di(D − µiI)−1z.


La dimostrazione del teorema puo essere trovata in [21]. Accenniamo quisoltanto a come nasce l’equazionef(µ) = 0, nota comeequazione secolare.Supponiamo per semplicita che siaρ = 1 e consideriamo la relazione

(D + zzT )q = µq,

doveµ rappresenta un autovalore diD+ ρzzT e q un corrispondente autovettore.Da questa segue(D − µI)q + z(zT q) = 0, da cuiq + (D − µI)−1z(zT q) = 0,e quindi

zT q + zT (D − µI)−1z(zT q) = 0.

In conclusione siamo arrivati a(1 + zT (D − µI)−1z)zT q = 0, ossia

f(µ)zT q = 0.

D’altra parte non puo esserezT q = 0; se infatti fosse cosı,q sarebbe autovettoredi D e quindi (essendoD diagonale) avrebbe un solo elemento diverso da 0 e diconseguenzaz dovrebbe avere un elemento uguale a 0. Si conclude pertanto chegli autovalori sono soluzioni dif(µ) = 0.

Per quanto riguarda le ipotesi del teorema, che a prima vistapossono sem-brare restrittive, possiamo fare alcune considerazioni. Il fatto chez non abbiacomponenti nulle e un’ipotesi generale. Sia infattizj = 0 per un indicej; alloradj e un autovalore di(D+ρzzT ) e il corrispondente autovettore e parallelo aej. Aquesto punto, si puo “sgonfiare” il problema, riconducendosi a uno di dimensionen − 1. Fra l’altro, come si osserva in [40], nei problemi pratici capita spesso cheil vettore z abbia alcune componenti nulle o molto piccole, e questo permette diaccelerare i calcoli proprio ricorrendo a tecniche di “sgonfiamento”. Questa con-siderazione ci porta naturalmente alla seconda ipotesi delteorema, quella secondocui D deve avere tutti elementi distinti. Cosa succede se questa condizione none verificata? L’idea per gestire questa situazione e dovuta a Dongarra e Sorensen[13] e viene applicata non solo sedj = dj+1 per qualche indicej, ma anche sela differenza|dj − dj+1| e piccola. In questa situazione si introduce forzatamenteuno 0 inz e si procede a sgonfiare il problema; essenzialmente, il procedimento eil seguente. Si considera una matrice di rotazioneGi ∈ R

2×2 tale che

GTi

(

zizi+1

)

=

(

τ0

)

,

ovviamente conτ =√

z2i + z2

i+1. Sia oraG = G(p, q, θ) la matrice di Givens in

cui (p, q) = (i, i + 1) eθ e quello diGi. E facile vedere che

GT (D + ρzzT )G = D + ρzzT + E,

dovez = (τ, 0)T , D = diag(d1, . . . , dn) con di = di e di+1 = di+1, eE e unamatrice tale che

‖E‖2 = |(di − di+1)cs|.

52 Capitolo 1

In particolare, sedi = di+1, D = D e la matriceE e la matrice nulla. Se|(di − di+1)cs| e minore di una soglia prefissata, si trascuraE e si sgonfia ilproblema lavorando suD e z. In [13] si dimostra che questo introduce un errorenel calcolo complessivo sufficientemente piccolo da non inficiare l’accuratezzadegli autovalori calcolati.

Si tratta a questo punto di vedere se e come l’equazione secolare puo essererisolta a poco costo e in modo stabile. Possiamo osservare a tale scopo che lafunzionef(µ) han poli (in d1, . . . , dn) ed e monotona fra i poli cosı che le radicisono alternate ai poli. All’inizio degli anni ’90 e stato individuato un procedi-mento iterativo per trovare le radici dif(µ) che, sfruttando queste caratteristichedella funzione, risulta convergere molto velocemente; unadiscussione esaurientesu questo punto si trova nell’articolo [3] del 1993.

Il costo computazionale del metodo Divide–and–conquer e analizzato in det-taglio in [40]. Qui ricordiamo soltanto che, grazie anche alla tecnica di sgon-fiamento a cui abbiamo accennato in precedenza, questo metodo e il metododa preferire quando si vogliano calcolare tutti gli autovalori e gli autovettori diuna matrice tridiagonale simmetrica. Se invece siamo interessati soltanto agliautovalori, allora puo convenire usare il metodo QR.

1.6.4 Metodo di bisezione

Questo metodo e particolarmente utile quando non si e interessati al calcolo ditutti gli autovalori ma soltanto ad alcuni di essi. Inoltre,il metodo non approssimagli autovettori, per i quali occorre eventualmente ricorrere a metodi diversi, comeil metodo di Iterazione inversa. Per introdurre il metodo dibisezione, ricordiamola definizione diinerziadi una matrice simmetrica:

Definizione 1.6.1.L’ inerzia di una matrice simmetricaA e la terna di numeriinteri In(A) = (n−, n0, n+), doven− e il numero di autovalori negativi,n0 ilnumero di autovalori nulli en+ il numero di autovalori positivi.

La base teorica del metodo e rappresentata dal seguente teorema, noto comeTeorema di inerzia di Sylvester [8]:

Teorema 1.6.3.SeA e una matrice simmetrica eX una matrice invertibile, alloraIn(A) = In(XTAX).

Dimostrazione. Sianon− e n′− il numero di autovalori negativi diA eXTAXrispettivamente, e siaN l’autospazio diA corrispondente ai suoi autovalori nega-tivi; allora dim(N) = n− e

xTAx < 0 per ogni x ∈ N.

Sia inoltreP l’autospazio diXTAX corrispondente aglin − n′− autovalori nonnegativi; allora dim(P ) = n− n′− e

xTXTAXx ≥ 0 per ogni x ∈ P.


EssendoX non singolare, anche il sottospazioXP ha dimensionen − n′− e laprecedente relazione puo essere riscritta come

xTAx ≥ 0 per ogni x ∈ XP.

Possiamo ora osservare che dim(N) + dim(XP ) = n− + n − n′−. Se fossen′− < n−, da questa relazione dedurremmo cheN ∩XP 6= ∅ e quindi esisterebbeun x ∈ R

n tale che contemporaneamentexTAx < 0 e xTAx ≥ 0. Essendoquesto ovviamente assurdo, deve esseren′− ≥ n−.

Rovesciando i ruoli diA eXTAX, con gli stessi ragionamenti si deduce chen− ≥ n′−, dimostrando cosı in conclusione chen− = n′−. Con analoghi passaggisi mostra che le due matrici devono avere lo stesso numero di autovalori positivi,e pertanto anche lo stesso numero di autovalori nulli.

Sulla base di questo teorema, possiamo ragionare cosı. Dato un qualunquenumero realec, supponiamo che esista una fattorizzazione

T − cI = LDLT ,

con D diagonale eL triangolare inferiore con elementi diagonali uguali a 1.Allora, in virtu del teorema 1.6.3 si ha

In(T − cI) = In(D).

Da questo segue che il numero di autovalori diT minori di c, uguale al numero diautovalori negativi diT − cI, e uguale al numero degli elementi diagonali diDminori di zero. La struttura diT rende molto semplice il calcolo di questo numero,che indicheremo d’ora in avanti comeσ(c). Infatti, essendoT tridiagonale, lamatriceL risulta bidiagonale e la relazioneT − cI = LDLT diventa

α1 − c β1

β1 α2 − c. . .

. . . . . . βn−1

βn−1 αn − c

=

1

l1. ... .. . . .

ln−1 1

d1. . .

.. .dn

1 l1. . . .. .

.. . ln−1

1

,

con ovvio significato dei simboli. Senza perdere in generalita possiamo supporrechec non appartenga allo spettro diT , e quindi che i coefficientid1, . . . , dn sianotutti diversi da zero. Sotto questa ipotesi possiamo calcolare questi coefficienti

54 Capitolo 1

nel modo seguente. Uguagliando gli elementi a destra e sinistra del segno di “=”,dalla prima riga si ricava

α1 − c = d1 e β1 = l1d1,

mentre dalle righe successive si ricava

αi − c = l2i−1di−1 + di e βi = dili.

Si ottiene allora la seguente relazione di ricorrenza:

{

d1 = α1 − c

dk = (αk − c) − β2k−1

dk−1k = 2, . . . , n.

(1.43)

Nell’articolo [9], nell’ambito della costruzione del software LAPACK, si di-mostra che questa relazione di ricorrenza e stabile, nel senso che i segni dei coeffi-cientidk (che e tutto quello che ci interessa) vengono calcolati bene a prescinderedagli errori di arrotondamento e dal fatto che non si applicanessuna strategia dipivoting. Viene inoltre svolta un’analisi approfondita della formula e degli ac-corgimenti che vanno utilizzati per una sua corretta applicazione, tenendo contodell’aritmetica floating point e della possibilita di overflow o underflow. In parti-colare, si considera l’eventualita che qualche coefficiente dk risulti uguale a zero,come sicuramente succede seT − cI non ammette la fattorizzazioneLDLT ;ebbene, nell’articolo citato si dimostra che sostituendo il coefficiente nullo con unnumero piccolo (positivo o negativo non ha importanza), si ottiene comunque unamatriceD con la stessa inerzia diT − cI. Vedremo piu avanti una spiegazione diquesto fatto.

Basandoci sulla relazione (1.43), possiamo calcolare il valore di σ(c) con ilseguente semplice algoritmo 1.6.1, il cui costo e di4n operazioni aritmetiche. Idati di cui l’algoritmo ha bisogno sono lo scalarec e la matriceT , espressa tramitela diagonaleα = (α1, . . . , αn) e la codiagonaleβ = (β1, . . . , βn−1).

Algoritmo 1.6.1. Algoritmo che calcola il numeroσ(c) di autovalori minori dicper una matrice tridiagonale simmetricaT .Dati: α, β, c;Risultato: s = σ(c);1. d = α1 − c;2. Sed < 0, alloras = 1, altrimentis = 03. Perk = 2, . . . , n

d = αk − c− β2k−1

dSed < 0, allora: s = s+ 1

Fine ciclo4. Fine


Questo algoritmo rappresenta uno strumento molto potente per calcolare gliautovalori, come vedremo fra poco. Ci piace osservare che esso si rivela moltoutile anche ad altri piu semplici scopi. Se, per esempio vogliamo sapere se lamatriceT e semidefinita positiva, basta calcolareσ(0): se risultaσ(0) = 0 dedu-ciamo che non ci sono autovalori diT minori di 0 e pertanto la risposta alla nostradomanda e affermativa.

L’algoritmo 1.6.1 e il cuore del metodo di bisezione, che puo prendere formediverse a seconda dello scopo specifico che ci si prefigge: calcolare un solo au-tovalore, calcolare un gruppo di autovalori consecutivi, calcolare gli autovaloricontenuti in un dato intervallo, ecc.

A titolo di esempio consideriamo il seguente problema: datoun intervallo[a, b] in cui cadono tutti gli autovalori diT ordinati in senso crescente daλ1 aλn,si vuol calcolareλj per unj fissato. L’intervallo[a, b] potrebbe essere

[a, b] = [−‖T‖1, ‖T‖1] = max1≤i≤n|αi| + |βi| + |βi−1| (1.44)

oppure (ricordando il Teorema 1.2.2 di Gerschgorin)

[a, b] = [min1≤i≤nαi − (|βi| + |βi−1|),max1≤i≤nαi + (|βi| + |βi−1|)], (1.45)

dove abbiamo postoβ0 = βn = 0. Se dividiamo in due l’intervallo[a, b], comefacciamo a sapere in quale dei due sottointervalli staλj? La risposta e semplice.Indichiamo conc il punto di mezzo dell’intervallo e calcoliamoσ(c): seσ(c) ≥ j,allora sappiamo cheλj e minore dic, altrimenti e maggiore o uguale dic. Suquesta base si puo impostare un algoritmo di bisezione comenell’algoritmo 1.6.2,dove i datikmax, ηa eηr rappresentano il numero massimo di iterazioni consentitee le tolleranze per il criterio di arresto.

Algoritmo 1.6.2. Algoritmo di bisezione per calcolare ilj-esimo autovalore apartire dal piu piccolo per una matrice tridiagonale simmetricaT .Dati: j, a, b, α, β, kmax, ηa, ηr;Risultati: λ ultima iterata calcolata

ind indicatore del risultatoind = 0: il criterio di arresto e stato soddisfattoind = 1: kmax iterazioni sono risultate insufficienti

1. Perk = 1, . . . , kmax

1. c = a+b2

2. calcolas = σ(c) tramite l’algoritmo 1.6.13. Ses ≥ j, allora: b = c, altrimenti:a = c4. Seb− a ≤ ηa + ηr|b|, allora:λ = c, ind = 0, Fermati

Fine ciclo2. ind = 13. λ = c4. Fine

56 Capitolo 1

Esempio 1.6.1.Consideriamo la matriceA definita da

ai,j = max(i, j) , i, j = 1, . . . , n.

Pern = 10, 20, . . . , 100, abbiamo trasformatoA in forma tridiagonale tramitel’algoritmo 1.4.1 e poi abbiamo calcolato l’autovaloreλ10 con l’algoritmo 1.6.2,usandoηa = ηr = 10−14. I risultati ottenuti sono in tabella 1.1, dovea e b sono

n λ10 a b k er

10 74.12158 −2.576e+01 1.032e+02 48 2.88e−1520 −5.006111 −1.143e+02 3.947e+02 55 3.33e−1530 −1.001643 −2.654e+02 8.683e+02 56 1.75e−1440 −1.710184 −4.789e+02 1.522e+03 57 5.06e−1550 −2.622955 −7.548e+02 2.355e+03 57 1.69e−1560 −3.739226 −1.093e+03 3.367e+03 57 6.89e−1570 −5.058758 −1.494e+03 4.556e+03 57 1.01e−1480 −6.581453 −1.957e+03 5.924e+03 57 4.72e−1590 −8.307265 −2.483e+03 7.468e+03 57 6.42e−15

100 −10.236170 −3.071e+03 9.189e+03 57 1.39e−15

Tabella 1.1 Esempio 1.6.1: metodo di bisezione

gli estremi dell’intervallo di Gerschgorin (1.45),k indica il numero di iterazionieffettuate eer e l’errore relativo fra il valore ottenuto con l’algoritmoe quellocalcolato con la funzioneeigdi Matlab. La tavola mostra che l’algoritmo funzionabenissimo per tutti i valori din considerati, come del resto ci si aspetta da unmetodo di bisezione. Piu precisamente, le ultime due colonne ci dicono che conun numero di iterazioni essenzialmente indipendente dan l’autovalore cercatoviene approssimato praticamente a precisione di macchina. �

Il metodo di bisezione nella formulazione che abbiamo qui presentato puoessere trovato nel libro di J.W. Demmel [8]. Ma nella maggiorparte dei libri dianalisi numerica, sia in inglese che in italiano, il metodo viene presentato in unaforma diversa. Piu precisamente gli autori seguono in generale la formulazioneoriginale proposta da J.W. Givens nel 1954 [20], nella qualesi deve supporre chela matriceT sia irriducibile. Questa formulazione fu poi studiata negli anni ‘60 daJ.H. Wilkinson ([4], [45]), che ne dimostro le ottime proprieta di stabilita, ma miseanche in luce alcuni problemi numerici che possono nascere nelle applicazioni.Per risolvere questi problemi, Wilkinson propose una modifica che essenzialmenteporta alla nuova versione. Raccontiamo questa storia nei dettagli, cominciandocon il descrivere il metodo nella formulazione di Givens.

Indichiamo conpk(λ) il polinomio caratteristico del minore principale di or-dine k di T per k = 1, . . . , n. Tenendo conto della struttura tridiagonale dellamatriceT −λI possiamo ricavare una relazione di ricorrenza fra questi polinomi.A questo scopo poniamo

p0(λ) ≡ 1.


Essendop1(λ) = α1 − λ, si vede banalmente che

p2(λ) = (α2 − λ)p1(λ) − β21p0(λ);

proseguendo con semplici calcoli si vede chep3 e legato ap2 ep1 dalla relazione

p3(λ) = (α3 − λ)p2(λ) − β22p1(λ),

analoga a quella che permette di esprimerep2 in termini di p1 e p0. Questa re-lazione e generalizzabile. Procedendo per induzione si vede infatti che vale laseguente relazione di ricorrenza:

p0(λ) ≡ 1p1(λ) = α1 − λpk+1(λ) = (αk+1 − λ)pk(λ) − β2

kpk−1(λ) , k = 1, . . . , n− 1.(1.46)

Il metodo di bisezione di Givens si basa su alcune proprietadella successione dipolinomi {p0, . . . , pn} che discendono dalla relazione (1.46) e sono indicate nelseguente teorema, noto cometeorema della successione di Sturm.

Teorema 1.6.4.Supponiamo che i coefficientiβ1, . . . , βn−1 siano tutti diversi dazero (T irriducibile). Allora per i polinomi definiti dalla relazione di ricorrenza(1.46) valgono le seguenti proprieta:

i) due successivi polinomipk epk−1 non hanno zeri in comune;

ii) gli zeri di pk−1 separano quelli dipk.

iii) Fissato un valorec, il numero di radici dipn minori di c e pari al numero divariazioni di segno nella successione numerica{p0(c), p1(c), . . . , pn(c)} 10.

Dimostrazione. Premettiamo che tutte le radici dei polinomip1, . . . , pn sonoreali, in quanto autovalori di matrici simmetriche, e quindi ha senso parlare diordinamento (e separazione).

i) La dimostrazione di questo punto e molto semplice. Se per qualche λ fossepk(λ) = pk−1(λ) = 0, in virtu della (1.46) avremmopk−2(λ) = pk−3(λ) =· · · = p0(λ) = 0, e questo e assurdo.

ii) Procediamo per induzione. La tesi e vera perk = 2. Infatti, l’unica radice dip1

e λ = α1 e le due radici dip2 sono date da(α1+α2)±

√(α1−α2)2+4β2

1

2 ; essendoβ1 6= 0, queste due radici sono una strettamente maggiore e l’altrastrettamenteminore diα1. Supponiamo ora che la tesi sia vera perk− 1. Osserviamo prima

10La presenza di uno 0 nella sequenza{p0(c), p1(c), . . . , pn(c)} deve essere contata come unae una sola variazione di segno. Ad esempio, nella sequenza{1,−1,−2, 0, 5} vanno contate 2variazioni di segno.

58 Capitolo 1

di tutto che la relazione (1.46) implica chepk e pk−2 hanno segno opposto inuno zero dipk−1, essendoβk−1 6= 0. Per ipotesi induttiva,pk−2 cambia segnouna e una sola volta in ogni intervallo delimitato da due successive radici dipk−1, e pertanto in ognuno di questi intervalli anchepk cambia segno un nu-mero dispari di volte. Di questi intervalli, ce ne sonok − 2; abbiamo quindiindividuato almenok − 2 radici di pk che separano quelle dipk−1. Le restantidue radici potrebbero in linea di principio trovarsi in qualcuno di questi inter-valli. Ma questo non e possibile. Infattipk e pk−2 hanno entrambi grado pario entrambi grado dispari, e si comportano quindi allo stessomodo a±∞. At-traversando la piu piccola radice dipk−1, pk e pk−2 hanno segno opposto, epertantopk deve cambiare segno un’altra volta per un valore diλ minore diquesta radice. Un ragionamento analogo ci fa vedere chepk ha una radice adestra della piu grande radice dipk−2. Questo conclude la dimostrazione delpuntoii).

iii) La dimostrazione di questo punto e molto tecnica e rimandiamo al libro [45] ilettori interessati ai dettagli.

Prima di tornare al metodo di bisezione vale la pena di evidenziare un’imme-diata conseguenza del teorema:

Corollario 1.6.1. Gli autovalori di una matrice tridiagonale simmetrica irriducibilesono tutti semplici.

Il punto iii) del precedente teorema fornisce l’idea di un algoritmo per cal-colareσ(c) diverso dall’algoritmo 1.6.1. Si tratta in questo caso di calcolare lasequenza

{1, p1(c), . . . , pn(c)}e contare le variazioni di segno in essa presenti.

Il problema insito in questa strada consiste nel fatto che, acausa della naturapolinomiale delle funzioni coinvolte, si puo incorrere troppo piu facilmente checon l’altra in fenomeni dioverflowo underflow. Questo fu chiaro appena si com-incio a studiare il metodo e Wilkinson e i suoi collaboratori misero in evidenza laquestione attraverso una ricca casistica nell’articolo [4] del 1967, dove si consid-era ad esempio la situazione in cui la matriceT abbia molti autovalori parecchiovicini fra loro. In questo caso il polinomiopn (che puo sempre essere fattorizzatonella formapn(λ) =

∏ni=1(λ−λi)) assume valori molto piccoli quando calcolato

in un puntoc vicino a questo gruppo di autovalori; da qui puo venire ununderflow.Nello stesso articolo si propone anche un rimedio a questo problema. Infatti nellasituazione precedente anchepn−1 ha un gruppo di radici molto vicine, perche lesue radici si alternano con quelle dipn; allora quandopn assume valori molto pic-coli, anchepn−1 fa altrettanto e insieme a lui anche un certo numero di polinomiprecedenti. Pertanto il valorepn(c)

pn−1(c)non e troppo piccolo. L’idea di Wilkinson e

dei suoi collaboratori fu quindi quella di usare le funzionirazionali

qk(λ) =pk(λ)

pk−1(λ)


invece dei polinomi per calcolare il numero di autovalori diT minori di un datoc. E evidente infatti che il numero di variazioni di segno in{p0(c), . . . , pn(c)}coincide con il numero di elementi negativi nella sequenza{q1(c), . . . , qn(c)}.

La “quadratura del cerchio” si ha ricavando dalla relazionedi ricorrenza(1.46) per i polinomipk una relazione per le funzioniqk; e facile vedere che,fissatoc, il valoreqk(c) e uguale adk per ognik = 1, . . . , n.

In conclusione, Wilkinson e i suoi collaboratori proponevano alla fine il meto-do nella versione che abbiamo visto all’inizio del paragrafo. Ma arrivavano a que-sta forma attraverso una strada molto piu complicata e inoltre, avendo alla base ilTeorema della successione di Sturm invece del teorema di inerzia di Sylvester, siportavano dietro l’ipotesi di irriducibilita che invece in linea teorica non serve.

Come gia abbiamo osservato parlando della relazione (1.43), il calcolo del-la sequenza{q1(c), . . . , qn(c)} non puo procedere seqj(c) e uguale a zero perun certo indicej, ovvero sepj(c) = 0. Il rimedio a questa situazione vieneofferto dal Teorema 1.6.4. Infatti, tenendo conto delle proprieta dei polinomip0, p1, . . . , pn, il numero di variazioni di segno nella sequenza{p0(c), . . . , pn(c)}non cambia se a un termine uguale a zero si sostituisce un numero ǫ diverso dazero, e quindi non cambia neppure il numero di termini negativi nella sequenza{q1(c), . . . , qn(c)}. Se i termini di questa sequenza vengono calcolati tramite la(1.43), occorre scegliereǫ abbastanza vicino a zero per avere la certezza che isegni dei termini successivi alj-esimo non vengano alterati.

Se si ha necessita di calcolare piu di un autovalore, si pu`o ovviamente ripeterel’algoritmo di bisezione tante volte quanti sono gli autovalori che ci interessano,ma questo modo di procedere e inutilmente dispendioso. In realta, l’algoritmo1.6.2 puo essere facilmente generalizzato in modo da calcolare un gruppo di auto-valori consecutivi, diciamo daλk aλj , conk ≤ j, senza ricominciare il processodi bisezione di[a, b] dall’inizio ogni volta.

Osserviamo infine che il metodo di bisezione e particolarmente adatto al-la realizzazione su elaboratori paralleli, come evidenziato in [21] e, piu diffusa-mente, nell’articolo [9].

1.6.5 Metodo di Jacobi

Il metodo di Jacobi e un metodo per calcolare tutti gli autovalori e gli autovet-tori di una matrice simmetrica. Esso si basa su principi molto semplici e risale aldiciannovesimo secolo [28], ma fu “recuperato” verso la finedegli anni ’50 delsecolo scorso perche facilmente realizzabile sugli elaboratori elettronici. Si inter-essarono alle proprieta di questo metodo molti analisti numerici, fra cui P. Henrici[24], J.H. Wilkinson [44] e H.P.M. van Kempen [41]. Nonostante oggi esistanometodi con velocita di convergenza molto piu elevata di quella del metodo di Jaco-bi, questo continua ad essere utilizzato per matrici di piccole dimensioni (nel qualcaso la velocita puo non costituire un grosso problema) a causa della sua sem-plicita che rende facile tradurlo in un algoritmo. Non solo, ma ci sono situazioniin cui questo metodo riesce a calcolare gli autovalori “piccoli” e i corrispondenti

60 Capitolo 1

autovettori con un’accuratezza molto maggiore dei metodi basati sulla tridiago-nalizzazione. Il motivo fondamentale di questo fatto e chein quelle situazioni latridiagonalizzazione introduce errori piu grandi di quanto non lo siano gli auto-valori. Questo aspetto e stato messo in luce in un articolo del 1992 di Demmel eVeselic [12] ed e alla base di un rinnovato interesse nel contesto del calcolo delladecomposizione ai valori singolari (SVD), di cui parleremonel capitolo successi-vo. Infine, un altro pregio di questo metodo, di cui parleremopiu avanti, e il suointrinseco “parallelismo”.

L’idea del metodo di Jacobi e quello di puntare alla diagonalizzazione dellamatrice simmetricaA mediante trasformazioni di similitudine ortogonali. Percapire i principi su cui queste trasformazioni si basano, supponiamo di volercalcolare autovalori e autovettori di una matrice2 × 2 simmetrica

A =

(

a ee b

)

, e 6= 0.

Ovviamente potremmo determinare gli autovalori risolvendo l’equazione caratte-ristica(a−λ)(b−λ)−e2 = 0 e poi, per quanto riguarda gli autovettori, calcolareuna soluzione non banale dei sistemiAx− λ1x = 0 eAx− λ2x = 0. Ma, cam-biando prospettiva, si puo effettuare il calcolo tutto in una volta. Si tratta di diago-nalizzare la matrice trovando una matrice ortogonaleQ tale che laB = QTAQ siadiagonale. Supponiamo di scegliere matrici di rotazione, ecerchiamo un valore diθ in modo tale che risultib21 = b12 = 0. Poiche risultab12 = ec2−es2+(a−b)cs,dobbiamo quindi risolvere l’equazione trigonometrica

e(cos2θ − sin2θ) + (a− b)cosθsinθ = 0,

eventualmente traformandola in modo da poter calcolare direttamentec e s senzapassare attraverso il calcolo diθ.

Per una matriceA di dimensione maggiore di 2, la diagonalizzazione nonpuo essere fatta in un solo passo e ci vorranno in generale infinite trasformazionidi similitudine per arrivare ad una matrice diagonale. Estendendo quanto appenadetto nel caso2×2, useremo matrici di Givens. In base a quanto gia osservato nelparagrafo 1.2.3, e facile vedere che, data una matrice di GivensG = G(p, q, θ),gli elementi della matriceB = GTAG coincidono con quelli diA sulle righe ecolonne di indici diversi dap e q, mentre gli elementi delle righep e q sono datidalle seguenti formule:

bpj = bjp = apjc− aqjsbqj = bjq = apjs+ aqjc

}

, j = 1, . . . , n, j 6= p, q

bpp = appc2 + aqqs

2 − 2apqscbqq = apps

2 + aqqc2 + 2apqsc

bpq = bqp = apq(c2 − s2) + (app − aqq)sc.

(1.47)

Supponiamo per il momento chep e q siano fissati in modo tale chep < q eapq 6= 0. Allora possiamo scegliereθ in modo che siabpq = 0, ovvero

apq(c2 − s2) + (app − aqq)sc = 0.


Ponendot = tgθ, dobbiamo allora risolvere l’equazione

apqt2 + (aqq − app)t− apq = 0, (1.48)

e poi calcolare

c =1√

1 + t2e s = tc. (1.49)

L’equazione (1.48) ha due soluzioni date da

t =−(aqq − app) ±

√

(aqq − app)2 + 4a2pq

2apq= −τ ±

√

1 + τ2,

conτ =

aqq − app

2apq.

Seτ = 0 le due soluzioni sono+1 e −1 (θ = ±π/4); altrimenti, esse sono invalore assoluto una minore di 1 e l’altra maggiore di 1. Come chiariremo piuavanti, nell’ambito del metodo di Jacobi occorre garantire|θ| < π/4; la sceltagiusta e pertanto

t =

{

−τ +√

1 + τ2 se τ ≥ 0

−τ −√

1 + τ2 altrimenti.

Per evitare errori di cancellazione conviene scrivere in modo diverso queste for-mule:

t =

{

1τ+

√1+τ2

se τ ≥ 01

τ−√

1+τ2altrimenti.

(1.50)

Chiediamoci ora quale effetto complessivo fa la trasformazione daA a B,al di la dell’annullarsi di un elemento extradiagonale. A questo scopo avremobisogno di utilizzare la norma di Frobenius delle matrici e alcune sue proprieta.Ricordiamo che questa norma e definita da

‖X‖F =

√

√

√

√

n∑

i,j=1

x2ij

per una generica matriceX ∈ Rn×n e che e invariante per trasformazioni di

similitudine. Per dimostrare questa ultima affermazione facciamo vedere prima ditutto che

‖X‖2F = tr(XTX).

PostoC = XTX, si ha

cjj =

n∑

i=1

xijxij =

n∑

i=1

x2ij, per j = 1, . . . , n;

62 Capitolo 1

da questo segue che

tr(C) =n∑

j=1

cjj =n∑

i,j=1

x2ij = ‖X‖2

F .

Sia oraY simile aX; alloraY eX hanno gli stessi autovalori e la stessa cosa valeperY TY eXTX. Poiche la traccia di una matrice e la somma degli autovalori, sideduce che tr(Y TY ) = tr(XTX), ossia‖Y ‖F = ‖X‖F .

Torniamo ora al nostro problema e indichiamo cone(A) e e(B) la sommadei quadrati degli elementi extradiagonali diA e di B, rispettivamente. Allorapossiamo scrivere

‖A‖2F =

n∑

i=1

a2ii + e(A), e ‖B‖2

F =

n∑

i=1

b2ii + e(B)

e dalla relazione di similitudine che intercorre fraA eB deduciamo che

e(A) − e(B) =

n∑

i=1

b2ii −n∑

i=1

a2ii. (1.51)

Entrambi i membri di questa uguaglianza possono essere semplificati. A sinistra,essendo invariati gli elementi delle righe e colonne con indici diversi dap e q, edessendoA eB simmetriche, resta soltanto

e(A) − e(B) = 2(a2pq − b2pq) + 2

n∑

j 6=p,q

[(a2pj + a2

qj) − (b2pj + b2qj)].

D’altra parte, dalle formule (1.47) e dal fatto chec2 + s2 = 1, segue

b2pj + b2qj = (apjc− aqjs)2 + (apjs− aqjc)

2 = a2pj + a2

qj ,

per ognij 6= p, q (ovvero gli elementi extradiagonaliapj e aqj cambiano ma lasomma dei loro quadrati resta invariata: se un elemento diB cresce rispetto alcorrispondente elemento diA, l’altro diminuisce di altrettanto). Essendo inoltrebpq = 0, otteniamo infine

e(A) − e(B) = 2a2pq. (1.52)

In altri termini, nel passaggio daA a B la somma dei quadrati degli elementiextradiagonali diminuisce di2a2

pq.Passiamo ora a esaminare il secondo membro della (1.51). Poiche peri 6= p, q

si haaii = bii, resta soltanto(b2pp + b2qq) − (a2pp + a2

qq). Otteniamo pertantol’uguaglianza

2a2pq = (b2pp + b2qq) − (a2

pp + a2qq).


Questo ci dice che quel che il decremento complessivo della parte extradiagonaleporta un equivalente incremento complessivo nella parte diagonale.

Torniamo alla relazione (1.52), da cui possiamo trarre uno spunto per la sceltadegli indici p e q. Infatti, nell’ottica di cercare la matriceG che “avvicini il piupossibile” la matriceB a una matrice diagonale, possiamo scegliere gli indicip eq in modo che sia

|apq| = maxi6=j|aij |.Questa scelta implica

a2pq ≥ e(A)

n(n− 1)(1.53)

e quindi, da (1.52),

e(B) ≤(

1 − 1

N

)

e(A), (1.54)

dove abbiamo postoN = n(n−1)2 .

Il metodo di Jacobi si basa sulle considerazioni fatte finora. Infatti, in questometodo si genera una successione di matrici{B(k)} con la regola

{

B(0) = AB(k) = GT

kB(k−1)Gk, k = 1, 2, . . .

doveGk = G(p, q, θ) e una matrice di Givens conp e q tali che

(b(k−1)pq )2 ≥ e(B(k−1))

n(n− 1),

ad esempio|b(k−1)

pq | = maxi6=j|b(k−1)ij |, (1.55)

eθ tale che|θ| ≤ π/4 e b(k)pq = b

(k)qp = 0. Con queste scelte, da (1.54) segue che

e(B(k)) ≤(

1 − 1

N

)k

e(A), (1.56)

ovvero la parte extradiagonale delle matriciB(k) tende a zero e sulla diagonaleemergono al crescere dik gli autovalori diA.

Per poter affermare che la successione{B(k)} converge, dobbiamo esseresicuri che esiste un ordinamento degli autovaloriλ1, . . . , λn tale che

limk→∞

b(k)ii = λi per i = 1, . . . , n.

Si puo dimostrare che questo ordinamento esiste se, come abbiamo supposto, adogni iterazione si sceglieθ tale che|θ| ≤ π/4; una traccia della dimostrazione sipuo trovare nel libro [32] di Parlett, mentre i dettagli si trovano in [19].

64 Capitolo 1

Una volta appurato che il metodo converge, possiamo chiederci quale e lavelocita di convergenza. La relazione (1.56) dice che la convergenza del metodoe lineare. Inoltre, poiche il fattore

(

1 − 1N

)

si avvicina a 1 al crescere din, ci siaspetta che la convergenza sia tanto piu lenta quanto piu grande e la dimensione diA. In realta, questa stima e abbastanza pessimistica, e si puo dimostrare [21] chela convergenza e asintoticamente quadratica sull’arco diN iterazioni, nel sensoche esiste una costanteγ > 0 per cui

e(B(k+N)) ≤ γe(B(k))

perk sufficientemente grande.Chiudiamo con un’ultima osservazione relativa al metodo diJacobi. Durante

un’iterazione un elemento extradiagonale viene annullatoe alcuni degli altri ven-gono modificati: in particolare, se un elemento interessatodalla trasformazionee nullo, il suo trasformato potrebbe essere (e in generale `e) diverso da zero.Per questo motivo la matriceA non viene tridiagonalizzata a priori: il metododistruggerebbe questa struttura, vanificando il lavoro fatto.

Il metodo di Jacobi e realizzato nell’algoritmo 1.6.3, chenecessita di alcuneprecisazioni e spiegazioni. Come al solito, nell’algoritmo si prevede per sicurezzaun massimo numero di iterazionikmax, oltre il quale il procedimento si arresta.Il criterio di arresto vero e proprio e basato sulla piccolezza die(B), doveB el’iterata corrente, relativamente alla parte extradiagonale originalee(A); dettaηla tolleranza, il criterio e pertanto il seguente:

e(B) ≤ ηe(A). (1.57)

Nell’algoritmo si evita di calcolare ogni voltae(B), ma si usa la relazione (1.52)per aggiornare ad ogni iterazione la grandezzaω, che rappresenta la meta die(B).L’aggiornamento viene fatto al passo 2.2, prima che la matriceB venga cambiata.

Osserviamo anche che nell’algoritmo usiamo gli algoritmi 1.2.2 e 1.2.3 pertrasformare la matriceB ad ogni iterazione invece delle formule (1.47). In-fine, non conserviamo le matriciUk = G1 . . . Gk; se si fosse interessati ancheal calcolo degli autovettori, si dovrebbe modificare l’algoritmo conservando leinformazioni essenziali per ricostruire queste matrici (cfr. [21]).

Algoritmo 1.6.3. Algoritmo di Jacobi per diagonalizzare una matrice simmetricaA ∈ R

n×n.Dati: n,A, η, kmax

Risultati: ind indicatore del risultatoind = 0: il criterio di arresto e stato soddisfattoind = 1: kmax iterazioni sono risultate insufficienti

k numero di iterazioni effettuateB ultima iterata calcolata


1.B = A, ω0 =∑

i<j a2ij

2. ω = ω0

2. Perk = 1, . . . , kmax

1. calcolap, q tali che|bpq| = maxi<j |bij|2. ω = ω − b2pq

3. τ =bqq−bpp

2bpq

4. Seτ ≥ 0, allora: t = 1τ+

√1+τ2

; altrimenti: t = 1τ−

√1+τ2

5. c = 1√1+t2

, s = ct

6. Perj = 1, . . . , nt1 = bpj , t2 = bqj

bpj = t1c− t2s , bqj = t1s+ t2cFine ciclo suj

7. Perj = 1, . . . , nt1 = bjp , t2 = bjqbjp = t1c− t2s , bjq = t1s+ t2c

Fine ciclo suj8. Seω < ηω0, allora: ind = 0, Fermati

Fine ciclo suk3. ind = 14. Fine

Il metodo che abbiamo finora descritto e chiamato metodo di Jacobiclassicoperche corrisponde alla forma originale del metodo. Si pu`o d’altra parte osser-vare che la scelta degli indicip e q secondo la (1.55) richiede i confronti fra tuttigli elementi extradiagonali allo scopo di individuare il piu grande in valore as-soluto. Anche tenendo conto della simmetria diA, il costo di questi confronti edell’ordine din2, decisamente inaccettabile, considerando anche che il costo del-la trasformazione effettiva e dell’ordine din. Per ovviare a questo inconvenientee stato proposto il metodo di Jacobiciclico nel quale gli elementi extradiagonalivengono annullati in un ordine prestabilito. Piu precisamente, supponendo di la-vorare nel triangolo superiore delle matrici, gli elementivengono annullati nelseguente ordine:

(1, 2) (1, 3) . . . (1, n)(2, 3) . . . (2, n)

...(n− 1, n)

Dopo un ciclo completo, composto daN iterazioni, si ricomincia. Ovviamente larotazione corrispondente a una coppia(p, q) viene saltata se l’elementobpq dellamatriceB corrente e nullo. Anche il metodo ciclico converge quadraticamente suun ciclo, come il metodo classico.

L’algoritmo 1.6.4 realizza il metodo ciclico. Esso e moltosimile all’algorit-mo 1.6.3, salvo il fatto che il numero massimo di iterazionikmax e sostituito daun numero massimo di ciclicmax.

66 Capitolo 1

Algoritmo 1.6.4. Algoritmo di Jacobi ciclico per diagonalizzare una matricesim-metricaA ∈ R

n×n.Dati: n,A, η, cmax

Risultati: ind indicatore del risultatoind = 0: il criterio di arresto e stato soddisfattoind = 1: cmax cicli sono risultati insufficienti

k numero di iterazioni effettuateB ultima iterata calcolata

1.B = A, k = 0, ω0 =∑

i<j a2ij

2. ω = ω0

3. Peri = 1, . . . , cmax

1. Perp = 1, . . . , n − 11. Perq = p+ 1, . . . , n

1. Se|bpq| >√

ωN

, allora:1. ω = ω − b2pq

2. τ =bqq−bpp

2bpq

3. k = k + 14. Seτ ≥ 0, allora: t = 1

τ+√

1+τ2; altrimenti: t = 1

τ−√

1+τ2

5. c = 1√1+t2

, s = ct

6. Perj = 1, . . . , nt1 = bpj , t2 = bqj

bpj = t1c− t2s , bqj = t1s+ t2cFine ciclo suj

7. Perj = 1, . . . , nt1 = bjp , t2 = bjqbjp = t1c− t2s , bjq = t1s+ t2cFine ciclo suj

8. Seω < ηω0, allora: ind = 0, FermatiFine scelta.

Fine ciclo suqFine ciclo sup

Fine ciclo sui4. ind = 15. Fine

Abbiamo inoltre inserito nell’algoritmo un controllo per saltare la rotazionese |bpq|, pur diverso da zero, risultasse piccolo; il criterio adottato e quello co-munemente accreditato in letteratura, ovvero

|bpq| ≤√

ω

N. (1.58)

In questo modo infatti la trasformazione viene effettuata soltanto se

b2pq >e(B)

n(n− 1),


da cui segue che la parte extradiagonale della matrice vieneridotta effettivamentealmeno di un fattore

(

1 − 1N

)

.Come vedremo dagli esempi, con il controllo (1.58) la variante ciclica diven-

ta in pratica equivalente a quella classica, nel senso che ilnumero di iterazionieffettive richieste per soddisfare un dato criterio di arresto e in generale non moltopiu grande di quello richiesto dal metodo classico per soddisfare lo stesso criterio.Stando cosı le cose, il costo (in termini di tempo di esecuzione) delle iterazioni inpiu e ampiamente compensato dal tempo risparmiato per la scelta dip e q.

Esempio 1.6.2.Consideriamo la matrice di Hilbert definita da

ai,j =1

i+ j − 1, i, j = 1, . . . , n,

Sen = 4 si ha e(A) ≃ 1.05 e gli autovalori, calcolati con la funzioneeig diMatlab, sono:9.6702e−5, 6.7383e−3, 1.6914e−1, 1.5002.

Jacobi classico Jacobi ciclicok (p, q) bpq ek (p, q) bpq ek

1 (1,2) 5.000e−01 7.8e−01 (1,2) 5.000e+00 7.8e−012 (1,3) 4.119e−01 5.2e−01 (1,3) 4.119e−01 5.2e−013 (1,4) 3.517e−01 1.5e−01 (1,4) 3.517e−01 1.5e−014 (2,3) 5.976e−02 1.2e−01 (2,3) 5.976e−02 1.2e−015 (2,4) 7.470e−02 5.3e−02 (2,4) 7.470e−02 5.3e−026 (1,2) 2.927e−02 3.3e−02 (1,2) 2.927e−02 3.3e−027 (1,4) −1.851e−02 2.0e−02 (1,3) −1.340e−02 2.7e−028 (1,3) −1.344e−02 5.8e−03 (1,4) −1.854e−02 5.8e−039 (3,4) 3.242e−03 3.5e−03 (2,3) 2.451e−03 4.6e−0310 (2,4) 1.849e−03 2.3e−03 (3,4) 3.236e−03 6.5e−0411 (1,3) 1.659e−03 6.2e−05 (2,3) −2.831e−04 5.1e−0412 (1,2) −2.711e−05 4.9e−05 (2,4) 3.558e−04 6.8e−0513 (1,4) −2.241e−05 3.7e−05 (1,2) −2.757e−05 5.6e−0514 (3,4) −1.890e−05 2.6e−05 (1,3) −3.118e−05 3.5e−0515 (1,3) 1.845e−05 2.7e−07 (1,4) −2.471e−05 8.4e−0716 (2,4) −1.860e−07 7.4e−08 (3,4) 5.953e−07 1.2e−09

Tabella 1.2 Esempio 1.6.2: risultati per la matrice di Hilbert di dimensione n = 4

Nella tabella 1.2 riportiamo i risultati ottenuti con il metodo di Jacobi classicoe con quello ciclico lavorando in doppia precisione, con tolleranza per il criteriodi arrestoη = 10−14. Per ogni iterazionek mostriamo gli indici(p, q) corrispon-denti, il valore dibpq, e infine il valore die(B) per la nuova matrice, indicato conek. Vediamo dalla tavola che entrambi i metodi riescono a soddisfare il criteriodi arresto in 16 iterazioni. In particolare, in virtu del controllo (1.58) la varianteciclica salta molte rotazioni “inutili”; vediamo ad esempio che le prime 6 itera-zioni risultano identiche a quelle del metodo classico. Entrambi i metodi danno

68 Capitolo 1

come risultato una matriceB i cui elementi diagonali approssimano con almeno13 cifre esatte gli autovalori diA. �

Esempio 1.6.3.Consideriamo la matrice dell’esempio 1.6.1, definita da

aij = max(i, j) , i, j = 1, . . . , n.

n Jacobi classico Jacobi ciclico incremento% in k risparmio%8 85 93 9.4 33.316 369 379 2.7 40.124 813 920 13.2 54.932 1459 1654 13.4 61.536 1831 2124 16.0 64.940 2275 2519 10.7 65.660 4946 5764 16.5 75.280 8621 10001 16.0 79.9100 13163 15009 14.0 82.9

Tabella 1.3 Risultati per l’esempio 1.6.3

In tabella 1.3 riportiamo per diversi valori din il numero di iterazioni richiestedalle due varianti del metodo per soddisfare il criterio di arresto conη = 10−14.Le ultime due colonne indicano in percentuale l’incrementonel numero di itera-zioni e il risparmio nel numero di operazioni che si hanno passando dalla varianteclassica a quella ciclica. Omologando operazioni aritmetiche e operazioni logiche,e indicando conk il numero di iterazioni, il costo del metodo classico e valutabilecome(12n+N)k operazioni, mentre quello del metodo ciclico e dato soltanto da12nk di operazioni. �

Il metodo di Jacobi ciclico ha una proprieta che lo rende molto interessantedal punto di vista algoritmico. Infatti si possono smpre individuare nell’insiemedelle rotazioni che compongono un ciclo completo dei gruppidi rotazioni che noninterferiscono l’una con l’altra (ad esempio, sen = 8, tali sono le rotazioni (1,2),(3,4), (5,6), (7,8)) e che quindi potrebbero essere eseguite contemporaneamentesenza cambiare il risultato. Sulla base di questa osservazione si puo scrivere unalgoritmo di Jacobiciclico parallelo, come raccontato in dettaglio in [21].

2DECOMPOSIZIONE AI VALORI

SINGOLARI

2.1 Introduzione

Un ruolo centrale in algebra lineare (teorica o numerica) egiocato dalle decompo-sizioni, o fattorizzazioni, delle matrici, che rappresentano uno strumento potenteper l’analisi dei problemi e per la progettazione di algoritmi risolutivi. Basti pen-sare alla fattorizzazioneLU che sta dietro al metodo di eliminazione di Gauss perla risoluzione di sistemi lineari, alla fattorizzazioneQR utile nella risoluzione diproblemi ai minimi quadrati lineari di rango pieno, alle fattorizzazioniLDLT edi Cholesky, solo per citare quelle piu comunemente usate.Anche nel campo deiproblemi di autovalori, un ruolo fondamentale e ricopertodalla diagonalizzazioneA = XΛX−1 quando esiste, e ancor piu dalla forma di SchurA = QTQH .

In questo capitolo ci occuperemo di una nuova decomposizione di matricidotata di proprieta teoriche e numeriche che ne fanno uno strumento assoluta-mente fondamentale in moltissime applicazioni. Si tratta della decomposizione aivalori singolari (SVD) di cui G. W. Stewart ha raccontato la storia nell’articolo[36]. La nascita della SVD viene fatta risalire al 1873, ad opera di E. Beltrami,che la derivo per matrici quadrate non singolari; l’estensione a matrici complessefu opera di L. Autonne [2] nel 1913, mentre l’estensione allematrici rettangolarie la generalizzazione delle principali proprieta della decomposizione risalgono al1936-39, ad opera di C. Eckart e G. Young([15],[16]). Come abbiamo gia fattonel precedente capitolo, anche qui ci occuperemo soltanto di matrici reali; inoltre,

Nel paragrafo 2.2 illustreremo la definizione e le principali proprieta dellaSVD per matrici inR

m×n e faremo un cenno a due importanti applicazioni: lacompressione di immagini e la risoluzione di problemi dei minimi quadrati linearirango-deficitari. Ci occuperemo poi nel paragrafo 2.3 dei metodi numerici per ilcalcolo della SVD. Dal momento che il problema e essenzialmente equivalente aopportuni problemi di autovalori per matrici simmetriche,i metodi che vedremosono varianti di quelli che abbiamo visto nel paragrafo 1.6.A questo proposito ri-cordiamo che prima di applicare un qualunque metodo per il calcolo di autovalorie/o autovettori di matrici simmetriche, ad eccezione del metodo di Jacobi, con-viene operare una trasformazione della matrice in forma tridiagonale superiore, aifini di rendere piu efficienti i procedimenti. Similmente, tutti i metodi per il cal-

70 Capitolo 2

colo della SVD, ad eccezione del metodo di Jacobi, prevedonoche la matrice siapreventivamente trasformata in una forma opportuna, che inquesto caso si rivelaessere la forma bidiagonale superiore, di cui ci occuperemonel paragrafo 2.3.2.

2.2 SVD e sue propriet a

2.2.1 Esistenza della SVD

Data una matriceA ∈ Rm×n e postop = min(m,n), una decomposizione ai

valori singolari (SVD) diA e una fattorizzazione della forma

A = UΣV T ,

doveU = (u1 · · · um) ∈ Rm×m e V = (v1 · · · vn) ∈ R

n×n sono ortogonalie Σ ∈ Rm×n e diagonale con elementi diagonaliσ1 ≥ σ2 ≥ . . . ≥ σp ≥ 0.Qui e nel seguito useremo il termine “diagonale” anche quandom 6= n nel sensoseguente: postoΣ = diag(σ1, σ2, . . . , σp), la matriceΣ ha la forma a blocchi

(

Σ0

)

se m > n e(

Σ 0)

se m < n.

Gli elementiσ1, . . . , σp sono chiamativalori singolari di A, le colonne diU sonoi vettori singolari sinistrie quelle diV sono ivettori singolari destri.

Osserviamo che per l’invarianza della norma−2 rispetto a trasformazioniortogonali, dalla SVD si deduce che

‖A‖2 = σ1. (2.1)

Da questo momento in poi supporremo semprem ≥ n; questa non e unarestrizione perche il caso opposto,m ≤ n, puo essere trattato come il primooperando con la matrice trasposta. La prima domanda da porcie se una SVD esisteper qualsiasi matrice. La risposta e affermativa, come si dimostra nel successivoteorema.

Teorema 2.2.1. Sia A ∈ Rm×n, con m ≥ n. Allora esistono due matrici

ortogonaliU ∈ Rm×m eV ∈ R

n×n tali che

UTAV =

(

Σ0

)

, (2.2)

conΣ = diag(σ1, σ2, . . . , σn) eσ1 ≥ σ2 ≥ . . . ≥ σn ≥ 0.

Dimostrazione. La dimostrazione procede per induzione. Pern = 1 em ≥ 1qualsiasi, la matriceA e di fatto un vettore inRm; in omaggio alla consuetudinedi indicare i vettori con lettere minuscole, chiamiamoa questo vettore. Poniamoσ1 = ‖a‖2 e indichiamo conH la matrice di Householder di dimensionem tale

DECOMPOSIZIONE AI VALORI SINGOLARI 71

cheHa = σ1e1. E evidente che abbiamo appena costruito una SVD diA, doveU = HT , Σ = σ1e1 eV = 1.

Sia oran > 1. Supponiamo che la tesi sia vera per matrici inRm×(n−1)

conm ≥ n − 1 qualsiasi e facciamo vedere che e vera per qualsiasi matriceA ∈ R

m×n. Per alleggerire la presentazione, nel seguito della dimostrazioneusiamo il simbolo‖ · ‖ senza pedici per indicare la norma−2 di vettori e matrici.

Consideriamo due vettoriu ∈ Rm e v ∈ R

n tali che‖u‖ = ‖v‖ = 1 eAv = ‖A‖u. Due siffatti vettori sicuramente esistono: possiamo prendere adesempiov unitario in modo che sia‖Av‖ = ‖A‖ (tale v esiste per definizionedi norme matriciali) e poiu = Av

‖Av‖ , da cui risulta immediatamente‖u‖ = 1

e Av = ‖Av‖u = ‖A‖u. Siano oraU ∈ Rm×m e V ∈ R

n×n due matriciortogonali tali che la prima colonna diU coincide conu e quella diV coincideconv, e poniamo

A = UTAV . (2.3)

Poiche

Ae1 = UTAV e1 = UTAv = UT ‖A‖u = ‖A‖UTu = ‖A‖e1,

la matriceA ha la forma

A =

(

σ1 wT

0 B

)

doveσ1 = ‖A‖, w ∈ Rn−1 eB ∈ R

(m−1)×(n−1). Se ora consideriamo il vettore

y =

(

σ1

w

)

, otteniamo

Ay =

(

σ21 + wTwBw

)

=

(

‖y‖2

Bw

)

e quindi‖Ay‖2

‖y‖2=

‖y‖4 + ‖Bw‖2

‖y‖2≥ ‖y‖2,

ovvero‖A‖2 ≥ ‖y‖2 = σ2

1 + wTw = ‖A‖2 + wTw.

Ma ‖A‖ = ‖A‖ percheA e ottenuta daA mediante trasformazioni ortogonali.Pertanto deve esserew = 0 e

A =

(

σ1 00 B

)

.

Per ipotesi induttiva la matrice B ammette una SVD, ovvero esistono due matriciortogonaliU1 ∈ R

(m−1)×(m−1) eV1 ∈ R(n−1)×(n−1) tali che

UT1 BV1 =

(

Σ2

0

)

72 Capitolo 2

con Σ2 = diag(σ2, . . . , σn) e ‖B‖ = σ2 ≥ . . . ≥ σn ≥ 0. Inoltre, da‖A‖ =‖A‖ = σ1 segue

σ21 = ‖A‖2 = ρ(AT A) = max(σ2

1 , ρ(BTB)) ≥ ρ(BTB) = ‖B‖2 = σ2

2

e pertantoσ2 ≤ σ1. Se allora definiamo

U =

(

1 00 U1

)

, V =

(

1 00 V1

)

e Σ =

(

σ1 00 Σ2

0 0

)

,

si ottiene la decomposizioneA = UΣV T

che, combinata con la (2.3), ci da

A = UΣV T , con U = U U e V = V V .

Resta cosı dimostrata l’esistenza di una SVD diA.

La SVD ha un’interpretazione geometrica molto semplice e interessante. In-dichiamo coneni e emi l’ i-esimo vettore della base canonica diR

n e Rm rispetti-

vamente, e osserviamo che

Avi = UΣV T vi = UΣeni = σiui (2.4)

per ogni i = 1, . . . , n. Osserviamo anche che le colonne diU e quelle diVcostituiscono una base ortonormale diR

m e di Rn rispettivamente. La relazione(2.4) ci dice allora che per ogni trasformazione lineare daR

n in Rm e possibile

individuare delle basi ortonormali dei due spazi tali che l’i-esimo vettore dellabase diRn viene trasformato in un multiplo dell’i-esimo vettore della base diR

m. Rispetto a queste basi, la trasformazione e identificata dalla matriceΣ, che abuona ragione puo essere vista come una matrice di strechting. Mediante la (2.4)si riesce a far vedere che ogni trasformazione lineare daR

n in Rm porta la sfera

unitaria diRn in un iperellissoide (i dettagli possono essere trovati in [40]).

Torniamo alla definizione di SVD. Se supponiamom > n, la matriceU puoessere partizionata nella formaU = (U1 U2) conU1 ∈ R

m×n eU2 ∈ Rm×(m−n),

da cui segue immediatamente la fattorizzazione

A = U1ΣV, (2.5)

talvolta chiamataSVD ridotta (in inglese,thin SVDoppuresingular value fac-torization ([37], [38])). La (2.5) ci dice che le ultimem − n colonne diU sono“ininfluenti” ai fini della SVD; piu avanti vedremo quale e il loro significato dalpunto di vista dell’algebra lineare.


2.2.2 SVD e decomposizioni spettrali

SeA e una matrice quadrata simmetrica esiste uno stretto legame fra valori sin-golari e autovettori e fra vettori singolari e autovettori.Infatti, essendoA ortogo-nalmente diagonalizzabile, possiamo scrivere

A = QΛQT = Q|Λ|sign(Λ)QT ,

con ovvio significato dei simboli. Poiche possiamo sempre supporre che gli ele-menti di |Λ| siano ordinati in senso decrescente, abbiamo di fatto una SVD di A,dove

U = Q, Σ = |Λ| e V T = sign(Λ)QT .

In altri termini i valori singolari sono i valori assoluti degli autovalori e i vettorisingolari sono autovettori.

SeA non e simmetrica i valori singolari e i vettori singolari non hanno legamidiretti con i suoi autovalori e autovettori, ma sono invece strettamente collegati aquelli delle due matrici simmetricheATA ∈ R

n×n eAAT ∈ Rm×m. Infatti da

A = UΣV T si deduce che

ATA = V Σ2V T e AAT = U

(

Σ2 00 0

)

UT (2.6)

e quindi entrambe le matrici ammettonoσ21 , . . . , σ

2n fra i loro autovalori; inoltre

AAT ha altrim − n autovalori nulli. Da (2.6) si deduce anche che le colonnedi V sono autovettori diATA e quelle diU lo sono perAAT ; in particolare, ivettori singolari “non essenziali”un+1, . . . , um sono autovettori corrispondentiall’autovalore nullo.

Dalle relazioni ora viste possiamo dedurre informazioni sull’unicita dellaSVD. Come prima cosa possiamo affermare che i valori singolari sono univo-camente determinati (in quanto autovalori di una matrice) equindi, in virtu del-l’ordinamento imposto, la matriceΣ e unica.

A causa dell’unicita dei valori singolari, e a dispetto della non unicita deivettori singolari di cui discuteremo fra un attimo, si e soliti parlare dellaSVD diuna matrice, e non di unaSVD.

Per quanto riguarda i vettori singolari, dalla definizione stessa di SVD dedu-ciamo che essi sono sempre definiti a meno del segno: seui evi sono una coppiadi vettori singolari associati all’i-esimo valore singolare, anche i loro opposti−ui

e−vi lo sono. Per esprimere questo concetto in formule consideriamo la SVD ri-dotta (2.5) e fissiamo una qualsiasi matrice di faseS ∈ R

n×n; allora e immediatoverificare che le matriciU1S, Σ eV S definiscono ancora una SVD ridotta.

Volendo approfondire il discorso, dobbiamo distinguere diverse situazioni.Supponiamo per semplicita che i valori singolari siano tutti diversi da zero. Allorapossiamo dedurre le seguenti considerazioni:

74 Capitolo 2

a) I vettori singolari sinistri “non essenziali”un+1, . . . , um non sono univoca-mente definiti, ma lo e il sottospazio da essi generato; questo discende dal fat-to che essi sono autovettori diAAT corrispondenti all’autovalore nullo che sipresenta con molteplicitam− n.

b) Se un valore singolareσi e semplice, l’essenza della non unicita dei vettorisingolari sinistro e destro corrispondentiui e vi sta tutta nella possibilita dicambiare il segno; questo discende dall’unicita degli autovettori normalizzaticorrispondenti ad autovalori semplici. Tenendo conto di questo, se tutti i valorisingolari sono semplici si dice che la SVD ridotta eessenzialmenteunica.

c) Se c’e un valore singolare ripetuto, i vettori singolariad esso corrisponden-ti non sono univocamente definiti ma lo e il sottospazio da essi generato; lamotivazione e la stessa del punto a).

Per riassumere le situazioni descritte, supponiamo che sia

σ1 = σ2 = · · · = σk > σk+1 > · · · > σn > 0.

Scelte comunque due matrici ortogonaliP ∈ Rk×k eZ ∈ R

(m−n)×(m−n), defini-amo

Q =

(

P 00 In−k

)

e Q =

(

Q 00 Z

)

;

allora si vede facilmente cheUQ, Σ eV Q individuano una nuova SVD diA.

Utilizzando la relazione fra valori singolari diA e autovalori diATA eAAT

si possono dedurre informazioni anche sul condizionamentodei valori singolari.Ricordando quanto detto nel paragrafo 1.3 per le matrici simmetriche, possia-mo infatti concludere subito che i valori singolari sono bencondizionati in sensoassoluto. In particolare, si puo dimostrare il seguente teorema [38].

Teorema 2.2.2.SiaA ∈ Rm×n con valori singolariσ1 ≥ σ2 ≥ . . . ≥ σp ≥ 0.

Siano inoltreσ1 ≥ σ2 ≥ . . . ≥ σp ≥ 0 i valori singolari della matrice perturbataA = A+ E. Allora

|σi − σi| ≤ ‖E‖2, i = 1, . . . , p.

Procedendo come abbiamo fatto nel paragrafo 1.3 a propositodegli autovaloridi matrici simmetriche e sfruttando la (2.1) si deduce che quandoA = fl(A),allora

|σi − σi|σi

≤ nσ1

σiǫm

e quindi nel processo di memorizzazione della matrice si conservano con unminimo di accuratezza soltanto i valori singolari per cui

σi ≥ nσ1ǫm. (2.7)

Ci preme sottolineare che questa relazione non dipende dalla “piccolezza” diσi

in assoluto, ma dalla sua “piccolezza” rispetto al valore singolare piu grandeσ1.


2.2.3 SVD e approssimazioni di rango basso

La SVD rappresenta uno strumento teorico molto potente per lo studio delletrasformazioni lineari e in particolare dei sottospazi fondamentali associati a unamatriceA, qualiN(A), R(A), N(AT ) eR(AT ). Alcuni risultati in questo sensosono raccolti nel successivo teorema, nel quale si mostra fra l’altro che la SVDe una decomposizionerank revealing, nel senso che da essa si deduce quale e ilrango diA.

Teorema 2.2.3.SiaA ∈ Rm×n conm ≥ n eA = UΣV T la sua SVD. Sia inoltre

r ≤ n il numero di valori singolari positivi, ossia

σ1 ≥ . . . ≥ σr > σr+1 = . . . = σn = 0. (2.8)

Allora:

i) PostoUr = (u1 · · · ur), Vr = (v1 · · · vr) eΣr = diag(σ1, . . . , σr), si ha

A = UrΣrVTr =

r∑

i=1

σiuivTi . (2.9)

ii) N(A) = span{vr+1, · · · , vn}.iii) R(A) = span{u1, · · · , ur} e pertanto il rango diA e esattamente uguale ar.

Dimostrazione.

i) Questo punto e banale: dalla SVD segueA =∑n

i=1 σiuivTi e quindi la (2.9) in

virtu della (2.8).ii) Siax ∈ N(A), cioe tale cheUΣV Tx = 0. EssendoU non singolare, questo

equivale a direΣV Tx = 0, ovvero (per l’ipotesi (2.8))ΣrVTr x = 0 ovvero

V Tr x = 0. In altri termini,x ∈ N(A) se e soltanto se e ortogonale av1, . . . , vr

ovvero se e soltanto se appartiene a span{vr+1, . . . , vn}.iii) Siax un vettore non nullo diRn ey = Ax. Da (2.9) segue

y = UrΣrVTr x = Urz

conz ∈ Rr, ovveroy ∈ span{u1, . . . , ur}. Viceversa, siay = Urz per qualche

z ∈ Rr. Poiche la matriceΣrV

Tr ha rango pieno, esiste sicuramente unx

diverso da zero inRn tale chez = ΣrVTr x, e quindiy ∈ R(A).

Ricordando che i sottospaziR(AT ) eN(AT ) sono i complementi ortogonali diN(A) eR(A) rispettivamente, dai puntiii) e iii) del teorema segue che

R(AT ) = span{v1, . . . , vr} e N(AT ) = span{ur+1, . . . , um}.

76 Capitolo 2

Uno degli utilizzi piu importanti della SVD nelle applicazioni nasce dalla(2.9). Essa ci dice che una matrice di rangor puo essere vista come combi-nazione lineare dir opportune matrici di rango 1. La (2.9) non e l’unico modoper esprimereA come combinazione di matrici di rango basso; ad esempio ognimatriceA = (a1 . . . an) ∈ R

m×n e banalmente somma din matrici di rango 1della forma(0 . . . ai . . . 0). L’importanza della (2.9) rispetto a tutte le altre possi-bili combinazioni, e che lak-esima somma parziale in (2.9) rappresenta la matricepiu vicina adA fra tutte le matrici di rango≤ k se si sceglie opportunamente ilmodo di misurare la distanza. Questo fatto e stabilito nel seguente teorema.

Teorema 2.2.4.SiaA ∈ Rm×n una matrice di rangor > 1 eA = UΣV T la suaSVD. Per ognik ∈ {1, . . . , r − 1}, siaBk = {B ∈ R

m×n : rango(B) ≤ k} eAk =

∑kj=1 σjujv

Tj . Allora

minB∈B‖A−B‖2 = ‖A−Ak‖2 = σk+1.

Dimostrazione.Come prima cosa facciamo vedere che vale la seconda uguaglian-za.Ak appartiene aBk per costruzione, e inoltre

‖A−Ak‖2 = ‖r∑

j=k+1

σjujvTj ‖ = ‖U ΣV T ‖2 = ‖Σ‖2 = σk+1,

doveΣ ∈ Rm×n e diagonale con elementi diagonali0, . . . , 0, σk+1, . . . , σn.

Sia oraB una qualunque matrice inBk eW = span{v1, . . . , vk+1}. Essendo

dim(N(B)) + dim(W ) ≥ n− k + k + 1 = n+ 1,

deve essereN(B)∩W 6= ∅. Siay ∈ N(B)∩W ; alloraBy = 0 ey =∑k+1

j=1 αivi,ovvero

V T y =

(

α0

)

con α = (α1, . . . , αk+1)T .

Senza perdere generalita possiamo supporre‖y‖2 = 1. Usando la definizione dinorma matriciale indotta e l’invarianza della norma−2 rispetto a trasformazioniortogonali, si ottiene

‖A−B‖2 ≥ ‖(A−B)y‖2 = ‖Ay‖2 = ‖UΣV T y‖2

= ‖ΣV T y‖2 ≥ σk+1‖V T y‖2 = σk+1‖y‖2 = σk+1

e quindi‖A−B‖2 ≥ σk+1, con il che la tesi e dimostrata.

Il teorema 2.2.4 ci dice che tramite la SVD si puo calcolare la migliore ap-prossimazione diA nel senso dei minimi quadrati nell’insiemeBk per un qualsiasi


k da 1 ar− 1. La bonta dell’approssimazione dipende dalla grandezza del primovalore singolare trascurato, essendo

‖A−Ak‖2

‖A‖2=σk+1

σ1. (2.10)

La compressione di immaginiUn’importante applicazione del teorema appena visto si trova nel cosiddetto

problema della “compressione delle immagini”, che puo essere sinteticamentedescritto nel modo seguente.

Per memorizzare un’immagine su un supporto digitale, la si deve scomporrein piccoli quadratini e associare ad ogni quadratino un certo numero di infor-mazioni utili a ricostruire la parte di figura contenuta nel quadratino stesso. Questeinformazioni sono valori numerici in genere ottenuti come media di informazionirelative a tutto il quadratino e sono memorizzate in un cosiddetto “pixel”(acronimodi “picture element”). Per immagini in bianco e nero ogni pixel e un numero cherappresenta la sfumatura di grigio media del quadratino corrispondente, andan-do da 0 (nero) a 1 (bianco). Per immagini a colori, tipicamente ad ogni pixelcorrispondono tre numeri che indicano l’intensita (media) dei tre colori primari:rosso, verde e blu. Pertanto un’immagine in bianco e nero pu`o essere memorizzatain una matricem×n, dovemn e il numero complessivo di quadratini; per un’im-magine a colori occorrono invece3mn valori. Piu grandi sonom e n, miglioree la qualita della ricostruzione dell’immagine nel momento della visualizzazione,ovvero la sua fedelta all’immagine originale. Valori tipici di m en vanno da circa500 per immagini semplici a 1024 o piu per immagini complesse.

Comprimere l’immagine significa ridurre il numero di informazioni neces-sarie per ottenere una buona ricostruzione dell’immagine stessa. Per “buona ri-costruzione” intendiamo che almeno a prima vista, grazie alle limitazioni dell’oc-chio umano, l’immagine approssimata risulta “uguale” a quella originale. Ebbene,il teorema 2.2.4 ci da uno strumento matematico concettualmente molto sempliceper raggiungere lo scopo. Consideriamo per semplicita immagini in bianco e nero.Se al posto della matriceA ∈ R

m×n dei pixels si memorizza la sua migliore ap-prossimazioneAk per unk fissato, l’occupazione di memoria passa dam × nlocazioni a(m+n)× k, corrispondenti aik vettoriu1, . . . , uk in R

m e i k vettoriσ1v1, . . . , σkvk inRn: scegliendo opportunamentek si puo ottenere un (notevole)risparmio di memoria. Conciliare bonta della ricostruzione e risparmio di memo-ria non e sempre facile; in particolare, non esiste una ricetta per sceglierek inmodo ottimale.

Esempio 2.2.1.Consideriamo l’immagine di figura 2.1 in alto a sinistra. Si trat-ta di un’immagine molto semplice, in cui non esistono sfumature di grigio, maunicamente il bianco e il nero, che e completamente descritta da una matriceA ∈ R

21×55 i cui elementi sono 0 o 1. Le altre immagini riportate in figurasonoquelle ricostruite usandoA1,A5 eA10: le prime due sono decisamente brutte, mala terza e accettabile. Per capire di piu cosa succede possiamo esaminare i dati

78 Capitolo 2

immagine originale 21 x 55

10 20 30 40 50

5

10

15

20

approssimazione di rango 1

10 20 30 40 50

5

10

15

20


10 20 30 40 50

5

10

15

20


10 20 30 40 50

5

10

15

20

Figura 2.1 la compressione di CIAO

k σk accuratezza risparmio%1 1.2654e+01 5.3274e−01 93.45 3.2530e+00 2.1067e−01 67.17 2.1165e+00 1.3805e−01 53.910 1.0729e+00 7.3445e−02 34.212 5.6479e−01 3.0306e−02 21.013 3.8350e−01 1.0387e−16 14.514 1.3143e−15 5.5772e−17 7.919 5.8416e−33 0 −−Tabella 2.1 Esempio 2.2.1: la compressione di CIAO

riportati in tabella 2.1. Nella seconda colonna di questa tabella riportiamo alcunivalori singolari diA calcolati dalla funzionesvd di Matlab; tenendo conto chei valori singolari daσ14 in poi sono molto piccoli, siamo autorizzati a dire che,numericamente parlando, la matrice ha rango 13 (infatti questa e la risposta chesi ottiene dalla funzionerank di Matlab). Nella terza colonna riportiamo l’accu-ratezza (2.10) con cuiAk approssimaA e nella quarta il risparmio di memoriache si ottiene considerandoAk al posto diA (il simbolo −− indica che non c’erisparmio perche(m + n)k > mn). Vediamo che la distanza relativa resta ab-bastanza grande perk ≤ 12 e diventa praticamente nulla (confrontabile con laprecisione di macchina) perk = 13, quandoAk praticamente coincide conA.Vediamo anche che il gioco vale la candela in termini di risparmio di memoria


soltanto sek ≤ 7. �

Esempio 2.2.2.Consideriamo l’immagine di figura 2.2 in alto a sinistra; essa edescritta da una matrice quadrata di dimensione 512 i cui valori singolari vannodaσ1 = 5.9894 × 104 aσ512 = 3.7542 × 10−2 diminuendo in modo abbastanzaregolare senza salti. Il rango della matrice calcolato con la funzionerank e 512.

Nella figura riportiamo le immagini ricostruite usandoA10, A100 e A250: la

immagine originale 512 x 512

100 200 300 400 500

100

200

300

400

500


100 200 300 400 500

100

200

300

400

500


100 200 300 400 500

100

200

300

400

500


100 200 300 400 500

100

200

300

400

500

Figura 2.2 la compressione di PUMPKINS

k σk accuratezza risparmio%1 5.9894e+04 2.0202e−01 99.610 3.9657e+03 6.2969e−02 96.1100 5.3487e+02 8.8083e−03 60.9200 1.6265e+02 2.6823e−03 21.9250 8.1367e+01 1.3318e−03 −−300 4.6314e+01 7.6673e−04 −−500 1.2165e+00 1.8031e−05 −−Tabella 2.2 Esempio 2.2.2: la compressione di PUMPKINS

prima e inaccettabile, mentre la seconda comincia ad essere accettabile anche se icontorni non sono ancora ben nitidi; la terza e decisamentemigliore, ma in questocaso non si ha nessun risparmio di memoria, come vediamo dalla tabella 2.2 chee organizzata come la tabella dell’esempio precedente. �

80 Capitolo 2

2.2.4 SVD, numero di condizionamento e pseudoinversa

SeA e una matrice quadrata non singolare, il suo numero di condizionamento inuna fissata norma naturale ek(A) = ‖A‖‖A−1‖. Se si utilizza la norma−2, dallarelazione (2.1) otteniamo un’espressione dik2(A) in termini dei valori singolari.Essendo infatti‖A‖2 = σ1 ed essendo i valori singolari diA−1 uguali ai reciprocidei valori singolari diA, si ha

k2(A) =σ1

σn. (2.11)

Vogliamo ora estendere la definizione dik2(A) a matrici non quadrate e/o rango–deficitarie. Per poter fare questo apriamo una parentesi perdimostrare un classicoteorema che da un’interpretazione importante del numero di condizionamento.Questo teorema afferma infatti che, qualunque sia la norma naturale scelta, ilreciproco dik(A) rappresenta la distanza relativa diA dall’insieme delle matricisingolari, misurata in quella norma. La dimostrazione del teorema segue da unaltrettanto importante lemma, noto come Lemma di Banach.

Lemma 2.2.1. Sia‖ · ‖ una norma naturale eM una matrice tale che‖M‖ < 1.Allora I −M e invertibile e

‖(I −M)−1‖ ≤ 1

1 − ‖M‖ . (2.12)

Dimostrazione.Dalla identita

(I −M)(I +M +M2 + · · · +Mk) = I −Mk+1

segueMk+1 = I − (I −M)Bk,

dove abbiamo postoBk = (I +M +M2 + · · · +Mk). D’altra parte, essendo‖M‖ < 1 si ha che‖Mk+1‖ ≤ ‖M‖k+1 tende a zero perk tendente all’∞. Seguepertanto che la successione{Bk} ha un limiteB tale cheI − (I −M)B = 0,ovveroB = (I −M)−1. Resta cosı dimostrata la prima parte della tesi. Per laseconda parte, basta osservare che le proprieta delle norme naturali implicano

‖Bk‖ ≤ ‖I‖ + ‖M‖ + ‖M‖2 + · · · + ‖M‖k =1 − ‖M‖k+1

1 − ‖M‖ ≤ 1

1 − ‖M‖per ognik. Passando al limite si ottiene la (2.12).

Siamo ora pronti a dimostrare il teorema annunciato sulla caratterizzazione delnumero di condizionamento.

Teorema 2.2.5.SiaA ∈ Rn×n non singolare e‖ · ‖ una norma naturale. Allora,

dettoB = {B ∈ Rn×n : detB = 0}, si ha

1

k(A)= minB∈B

‖A−B‖‖A‖ .


Dimostrazione.Per dimostrare la tesi e ovviamente sufficiente dimostrareche

min{‖δA‖ : A+ δA ∈ B} =1

‖A−1‖ .

SiaδA tale che‖δA‖ < 1‖A−1‖ . Allora ‖A−1δA‖ ≤ ‖A−1‖‖δA‖ < 1 e il Lemma

di Banach implica cheI+A−1δA e non singolare, ovveroA+δA e non singolare.Pertanto

min{‖δA‖ : A+ δA ∈ B} ≥ 1

‖A−1‖ .

Costruiamo ora una matriceδA tale che‖δA‖ = 1‖A−1‖ e A + δA e singo-

lare. Dalla definizione di norma naturale, esiste unx ∈ Rn tale che‖x‖ = 1

e ‖A−1‖ = ‖A−1x‖. Poniamoy = A−1x‖A−1x‖ = A−1x

‖A−1‖ , da cui segue‖y‖ = 1, e

definiamoδA = −xyT

‖A−1‖ . Allora si ha

‖δA‖ = maxz 6=0‖xyT z‖

‖A−1‖‖z‖ =‖x‖

‖A−1‖maxz 6=0|yT z|‖z‖

e, applicando la disuguaglianza di Cauchy-Schwartz,‖δA‖ = 1‖A−1‖‖y‖ = 1

‖A−1‖ .Inoltre,A+ δA e singolare perche per definizione diy si ha

(A+ δA)y = Ay − xyT y

‖A−1‖ =x

‖A−1‖ − x‖y‖‖A−1‖ = 0.

Alla luce di questo teorema viene naturale estendere la definizione di numero dicondizionamento a matrici rettangolariA ∈ R

m×n di rango pieno. Possiamoinfatti definirek(A) comeil reciproco della distanza relativa diA dall’insiemedelle matrici inR

m×n rango-deficitarie. Se specifichiamo la norma scegliendo lanorma–2, il teorema 2.2.4 ci permette di dare una forma ak2(A):

k2(A) =σ1

σp,

conp = min(m,n). Piu in generale, ser > 1 e il rango diA, possiamo definirenello stesso spirito

k2(A) =σ1

σr, (2.13)

il cui inverso rappresenta la distanza relativa diA dall’insieme delle matrici inR

m×n di rango inferiore ar. Infine, seA e una matrice di rango 1, possiamoancora definirek(A) come inverso della distanza relativa diA dalla matrice nulla,ovverok(A) = 1, qualunque norma si utilizzi.

Il numero di condizionamentok2(A) definito da (2.13) per matrici di rangomaggiore di 1 e da 1 per matrici di rango 0 ha anche un’altra interpretazione intermini della matrice pseudoinversa di Moore-Penrose che andiamo a definire.

82 Capitolo 2

Definizione 2.2.1.SiaA ∈ Rm×n di rangor eA = UΣV T la sua SVD. Indichi-

amo conΣ+ ∈ Rn×m la matrice diagonale con elementi diagonali uguali a1

σiper

i = 1, . . . , r e0 peri = r + 1, . . . ,min(m,n). Allora, la matrice

A+ = V Σ+UT

e dettapseudoinversa di Moore-Penrosedi A.

L’appellativo “pseudoinversa” nasce dal fatto, facilmente verificabile, che seA e quadrata e non singolare, alloraA+ = A−1. Inoltre, qualunque siaA, lamatriceA+ e l’unica matrice che soddisfa le seguenti proprieta:

(a) AA+A = A(b) A+AA+ = A+

(c) (AA+)T = AA+

(d) (A+A)T = A+A.

E assolutamente banale verificare che

k2(A) = ‖A‖2‖A+‖2

qualunque siano le dimensioni e il rango diA.

2.2.5 Rango numerico

Dal teorema 2.2.4 discende immediatamente il seguente corollario:

Corollario 2.2.1. SiaA ∈ Rm×n conm ≥ n di rango pieno. Allora, ogni matrice

B ∈ Rm×n tale che‖A−B‖2 < σp ha rango pieno.

In altre parole, l’insieme della matrice di rango pieno e aperto e denso inR

m×n. Questo implica che, introducendo piccole perturbazioni in una matrice dirangor < p si ottiene quasi certamente una matrice di rango pieno. Possiamoquindi concludere che in aritmetica floating point quasi tutte le matrici hanno tuttii valori singolari maggiori di zero, e questo fa diventare problematico il calco-lo del rango di una matrice, o semplicemente stabilire se ha rango pieno o no,basandosi sul teorema 2.2.3. Questo non ci sorprende: e la stessa situazione chesi presenta quando si vuol decidere se una matrice quadrata `e singolare o no. Inquesta situazione ci si deve accontentare di una risposta insenso numerico: unamatrice enumericamente singolarese il suo numero di condizionamentok(A)in una qualche norma e maggiore di una soglia prefissata (di solito pari o pro-porzionale all’inverso della precisione di macchinaǫm), ovvero se la sua distanzadall’insieme delle matrici singolari e molto piccola.

Allo stesso modo diremo che una matrice rettangolare enumericamente rango-deficitariaseσn e minore di una certa soglia dipendente daǫm.



A =

1 2 312

13

56−1 4 3

che ha rango 2, essendo la terza colonna somma delle prime due. Definendola matrice in Matlab e calcolandone i valori singolari tramite la funzionesvdsiottengono i valori6.1420, 1.8252 e 3.8639e−16. Il rango della matricefl(A)e 3, ma il terzo valore singolare e chiaramente trascurabile e si puo dire che,numericamente parlando, la matrice e rango deficitaria. �

Piu in generale, possiamo introdurre il concetto dirango numericodi unamatrice. Fissata una tolleranzaǫ dipendente daǫm, diremo rango numerico diuna matrice il numero dei suoi valori singolari maggiori diǫ. Questa definizioneda informazioni precise quando c’e un salto netto fra i valori singolari maggioridi ǫ e quelli piu piccoli, come succede nell’esempio precedente o nell’esempio2.2.1 se prendiamo una qualunque soglia “ragionevole”. La situazione e piu deli-cata quando questo salto non c’e, nel qual caso l’informazione puo essere ancoravalida qualitativamente ma non quantitativamente. Mostriamo questa situazionenell’esempio seguente.

Esempio 2.2.4.Definiamo la matriceA ∈ Rn×n diagonale, conaii = 0.5i per

i = 1, . . . , n. Essendoσi = aii per ognii, la matrice e non singolare qualunquesian. D’altra parte 1

k2(A) = 0.51−n diminuisce al crescere din; supponendo dilavorare in doppia precisione con locazioni di memoria di 64bits, la precisione dimacchina eǫm = 0.5−52 e pertantoA e da considerare numericamente singolarepern > 52. Cosa possiamo dire dal punto di vista del rango numerico? Esaminia-mo i valori singolari daσ41 aσ52 calcolati da Matlab (senza errori perche trattasidi potenze di 2):

σ41 = 4.5475e− 13, σ42 = 2.2737e− 13 σ43 = 1.1369e− 13σ44 = 5.6843e− 14, σ45 = 2.8422e− 14 σ46 = 1.4211e− 14σ47 = 7.1054e− 15, σ48 = 3.5527e− 15, σ49 = 1.7764e− 15σ50 = 8.8818e− 16, σ51 = 4.4409e− 16, σ52 = 2.2205e− 16.

Fissandoǫ = 10−14 avremo la risposta che il rango numerico diA e 46; ma perǫ = 5 × 10−14, la risposta sara diversa e il rango numerico sara uguale a44. �

In Matlab il rango numerico di una matrice puo essere calcolato tramite lafunzionerank, a cui abbiamo fatto riferimento anche negli esempi 2.2.1 e 2.2.2.Questa funzione usa la soglia

ǫ = max(m,n)σ1ǫm (2.14)

e per la matrice dell’esempio precedente calcola rango numerico uguale a 47 pern = 50, 44 pern = 500, 43 pern = 1000 e 42 pern = 2000.

84 Capitolo 2

La scelta (2.14) e assolutamente ragionevole se si pensa che i valori singolaripiu piccoli di nσ1ǫm sono quasi sicuramente affetti da grossi errori relativi (cfr.(2.7). Il fattore max(m,n) e una salvaguardia per tenere conto anche degli erroriintrodotti durante il calcolo dei valori singolari. Ma questa scelta e ragionevoleanche da un altro punto di vista. Sapendo che il rango di una matrice non cam-bia se la moltiplichiamo per uno scalare, possiamo pretendere che anche il rangonumerico non cambi, o cambi molto poco, nella stessa situazione. Ebbene, semoltiplichiamo la matriceA dell’esempio precedente per una costante, ad esem-pio per10 o 10−30 o 1030, la funzionerank da ancora le stesse risposte, con unleggero scarto pern = 50 (48 invece che 47 in un caso), a dispetto della grandezzaassoluta dei valori singolari non “contati”.

Abbiamo appena visto che la SVD e uno strumento potente per calcolare ilrango (numerico) di una matrice. Ad oggi e di fatto il piu affidabile strumento attoa questo scopo. D’altra parte, come vedremo, il calcolo della SVD e in generalemolto costoso, potremmo dire troppo costoso soprattutto per matrici di grandi di-mensioni con rangor molto minore di min(m,n). Questo ha spinto molti analistinumerici a cercare dei metodi alternativi basati su fattorizzazioni diverse, il cuicalcolo richieda un minor numero di operazioni. La prima idea in questo sensorisale agli anni ’60 del secolo scorso, ed e dovuta a Golub, il quale propose unaforma di fattorizzazione QR con pivoting che in molti casi permette di individ-uare il rango della matrice osservando la “coda” della matrice triangolareR (cfr.[21]). Ma questo algoritmo non e completamente affidabile (esiste un famosocontroesempio dovuto a Kahan, riportato anche sul libro stesso di Golub e VanLoan), ragion per cui molti autori hanno proposte altre rank-revealing fattoriz-zazioni, comunque sempre basate sulla fattorizzazione QR.Questo argomento, siapure affascinante, non rientra negli scopi di questi appunti; gli interessati possonoleggere, ad esempio, gli articoli [6] e [27].

2.2.6 SVD e problema dei minimi quadrati lineari

Uno dei problemi piu importanti dell’algebra lineare e ilcosiddetto problema diminimi quadrati lineari, nel quale, datiA ∈ R

m×n eb ∈ Rm, si deve calcolare un

vettorex ∈ Rn soluzione del problema

minx‖Ax− b‖2. (2.15)

Le principali proprieta del problema sono descritte e dimostrate, ad esempio, in[5] o [21]. Qui ricordiamo soltanto i fatti principali:

1) l’insieme delle soluzioniX e un insieme non vuoto, chiuso e convesso;2) la soluzione e unica se e soltanto ser = n;3) esiste un solo vettorex∗ ∈ X di norma euclidea minima, ovvero tale che‖x∗‖2 = minx∈X‖x‖2.

Il seguente teorema mette in relazione il problema (2.15) con la SVD dellamatriceA.


Teorema 2.2.6.SiaA ∈ Rm×n, conm > n, di rangor, e siaA = UΣV T la sua

SVD. Sia inoltreb ∈ Rm. Allora la soluzione di minima normax∗ del problema

(2.15)e data da:

x∗ =r∑

i=1

uTi b

σivi. (2.16)

Dimostrazione.Da‖Ax− b‖2 = ‖UΣV Tx− b‖2 = ‖ΣV Tx−UT b‖2, ponendoz = V Tx, si ottiene

‖Ax− b‖22 = ‖Σz − UT b‖2

2

=∑n

i=1(σizi − uTi b)

2

=∑r

i=1(σizi − uTi b)

2 +∑n

i=r+1(uTi b)

2 ≡ ψ(z).

Il funzionaleψ(z) e minimizzato ponendo

zi =uT

i b

σiper i = 1, . . . , r.

Se r = n, abbiamo in questo modo completamente definito l’unico punto diminimo diψ e quindi l’unica soluzione di (2.15) che sara data da (2.16).Ser < n, le componentizr+1, . . . , zn possono essere scelte in un modo qualsiasi,dando luogo a tutti i possibili punti di minimo diψ e quindi, tramite la relazionex = V z, a tutte le soluzioni di (2.15). Tenendo conto del fatto che‖x‖2 =‖V z‖2 = ‖z‖ per ogniz, la soluzione di minima norma di (2.15) si ottiene incorrispondenza del punto di minimo diψ di norma minima, che e ovviamente

z∗ = (uT

1 b

σ1, . . . ,

uTr b

σr, 0, . . . , 0)T .

A questo corrispondex∗ = V z∗, ovvero (2.16).

Grazie a questo teorema e immediato verificare chex∗ = A+b, doveA+ e lapseudoinversa di Moore-Penrose diA che abbiamo definito nel paragrafo 2.2.4.

Ci piace completare questo paragrafo con un cenno all’analisi del condiziona-mento del problema dei minimi quadrati lineari, in cui intervienek2(A). Esistonomolti risultati in questo senso (tutti riferiti al caso di rango pieno) che differisconoper particolari tecnici, ma non nella sostanza; riportiamoqui il seguente teoremadi perturbazione tratto da [8].

Teorema 2.2.7.SiaA ∈ Rm×n, conm ≥ n, una matrice di rango pieno, e sia

x∗ 6= 0 la soluzione del problema dei minimi quadrati lineari (2.15) con residuor∗ = Ax∗ − b. Siano inoltreδA ∈ R

m×n e δb ∈ Rm tali che

ǫ ≡ max(‖δA‖2

‖A‖2,‖δb‖2

‖b‖2) <

1

k2(A).

86 Capitolo 2

Allora, dettay∗ una soluzione del problema perturbato

miny‖(A+ δA)y − (b+ δb)‖2,

si ha che:

i) A+ δA ha rango pieno,ii) L’errore relativo soddisfa la disuguaglianza

‖x∗ − y∗‖2

‖x∗‖2≤ ǫ

(

2k2(A)

cosθ+ tanθ · k2(A)2

)

+O(ǫ2) ≡ ǫkLS +O(ǫ2) (2.17)

dove sinθ = ‖r∗‖2

‖b‖2.

La limitazione (2.17) prospetta tre diverse situazioni. Laprima e quandoθ e uguale a 0 o molto piccolo (ovvero il residuor∗ e nullo o molto piccolo);in questo caso il condizionamento della soluzione e essenzialmente misurato da2k2(A). La seconda situazione e quandoθ non e vicino a zero e neppure vicinoa π/2; in questo caso il condizionamento della soluzione e peggiore di prima,in quantokLS e dell’ordine di grandezza dik2(A)2. Infine, la terza situazione equandoθ e vicino aπ/2 (ovverox∗ e vicino a zero), nel qual casokLS diventagrandissimo, al limite infinito, anche sek2(A) e piccolo. Queste osservazioni,sia pur qualitative, ci dicono che la sensibilita della soluzione di un problema aiminimi quadrati lineari dipende non solo daA, come succede nel caso dei sistemilineari, ma anche dal vettoreb.

Osservazione 2.2.1.Vogliamo concludere questo paragrafo con un’osservazioneriguardante il software per la risoluzione dei problemi deiminimi quadrati linea-ri. In ambiente Matlab troviamo la funzionelinsolve, che non usa la SVD ma ilmetodo QR con pivoting [21] e segnala con unwarningil caso in cui la matrice siarango-deficitaria. Questa scelta e probabilmente motivata dal fatto che il calcolodella SVD e solitamente abbastanza costoso e, d’altra parte, il metodo QR conpivoting riesce a risolvere la maggior parte dei problemi. Se invece passiamo aLAPACK, troviamo diverse alternative: il driver xGELSY usauna fattorizzazioneortogonale completa rank-revealing che viene ottenuta a partire da una QR conpivoting; i drivers xGELSS e xGELSD usano invece la SVD, calcolata con duemetodi diversi.

2.3 Metodi numerici per il calcolo della SVD

2.3.1 Caratteristiche generali dei metodi

I risultati del paragrafo 2.2.2 sui legami fra i valori/vettori singolari diA e gliautovalori/autevettori diAAT eATA, indicano una possibile strada per calcolarela SVD di una matrice. I passi potrebbero essere i seguenti:


1. CalcolareC = ATA;2. DiagonalizzareC, ovvero calcolare la fattorizzazioneC = V ΛV T , con V

ortogonale eΛ = diag(λ1, . . . , λn) conλ1 ≥ λ2 ≥ . . . ≥ λn;3. CalcolareΣ = Λ

1

2 ;4. Calcolare le colonne diU dalla relazioneAV = UΣ, ovveroui = Avi

σiper

i = 1, . . . , r e poiur+1, . . . , um come completamento della base ortogonale diR

m.

Questa strada, a cui Stewart in [38] da il nome di algoritmocross-product,e poco costosa in confronto alle alternative che vedremo successivamente, maha grosse limitazioni dal punto di vista dell’accuratezza.Per capire questa af-fermazione, limitiamo l’analisi ai valori singolari. Supponiamo per semplicitadi non commettere errori nel calcolo diC e consideriamo soltanto l’effetto del-la sua memorizzazione in precisione finita. SiaC = fl(C) con autovaloriλ1,λ2, . . ., λn . Sappiamo cheλi puo perdere completamente accuratezza rispettoa λi seλi < nǫmλ1. Ma, ricordando cheλi = σ2

i , questo equivale a dire chel’approssimazione diσi per questa strada puo essere completamente inaccurata se

σi < n√ǫmσ1.

Questa limitazione e decisamente inaccettabile, soprattutto pensando che se nonsi passa daA aC si perde accuratezza nei valori singolari per cui

σi < nǫmσ1.

In vista di queste osservazioni, i metodi effettivamente usati per il calcolo del-la SVD sono varianti dei metodi per calcolare autovalori e autovettori di matricisimmetriche, nei quali si “mima” l’algoritmo cross-product evitando il calcolo es-plicito di C e si sfrutta la struttura della matrice per ottenere la massima efficienzapossibile. Di quale struttura stiamo parlando? Ricordiamoche tutti i metodi peril calcolo di autovalori/autovettori di matrici simmetriche richiedono la preventi-va trasformazione della matrice in forma tridiagonale. In questo caso la matricee ATA e lo scopo puo essere ottenuto trasformando preventivamente A in unamatrice bidiagonale superiore:

A→ B =

(

B0

)

, con B ∈ Rn×n bidiagonale superiore. (2.18)

Vedremo nel prossimo paragrafo che questa trasformazione puo essere ef-fettuata in modo stabile usando trasformazioni di Householder. Inoltre, essa nonpresenta svantaggi dal punto di vista della accuratezza perche vale il seguente teo-rema (cfr. [8, Teorema 5.13]) che sancisce il buon condizionamento dei valorisingolari di una matrice bidiagonale.

Teorema 2.3.1.SianoB e B due matrici bidiagonali, con elementi diagonaliα1, . . . , αn e α1, . . . , αn rispettivamente, e elementi sopradiagonaliβ1, . . . , βn−1

88 Capitolo 2

e β1, . . . , βn−1 rispettivamente, doveαi = αi(1 + ǫi) e βi = βi(1 + θi) con|ǫi|, |θi| ≤ η << 1. Siano inoltreσ1, . . . , σn i valori singolari diB e σ1, . . . , σn

quelli di B. Allora

i) da σi = 0 segueσi = 0;ii) da σi 6= 0 segue|σi−σi|

σi≤ (4n − 2)η + O(η2).

Cio detto, tutti i metodi per il calcolo della SVD, ad eccezione del metodo diJacobi, seguono lo schema seguente:

1. Calcolare una fattorizzazioneA = UBBVBT conB bidiagonale superiore e

UB eVB ortogonali (bidiagonalizzazione diA);2. Calcolare la SVD diB, ossiaB = UΣV T ;3. Combinare le due fattorizzazioni per ottenere la SVD desiderata (A = UΣV T

conU = UBU eV = V TB V ).

Il nucleo computazionale dei procedimenti diventa quindi il punto 2. nelquale si utilizzano metodi iterativi ad hoc per il calcolo della SVD di matricibidiagonali. Con riferimento alla (2.18) osserviamo a questo proposito che esufficiente lavorare sullaB. Infatti, seB = U ΣV T , la SVD di B e data daB = UΣV T , con

U =

(

U 00 Im−n

)

, Σ =

(

Σ0

)

, V = V .

2.3.2 Riduzione in forma bidiagonale

Qualsiasi matriceA ∈ Rm×n, conm ≥ n, puo essere trasformata in forma bidi-

agonale superiore tramite trasformazioni ortogonali. Infatti esistono sicuramentedue matriciUB ∈ R

m×m e VB ∈ Rn×n ortogonali e una matrice bidiagonale

superioreB ∈ Rm×n tali che

UTBAVB = B. (2.19)

La dimostrazione di questo fatto e costruttiva. Consideriamo infatti la sequenzadi trasformazioni

A → U1A→ U1AV1 → U2U1AV1 → U2U1AV1V2 → . . .

Un−2 · · ·U1AV1 · · · Vn−2 → Un−1 · · ·U1AV1 · · ·Vn−2 →

Un · · ·U1AV1 · · ·Vn−2 ≡ B,

doveUk, perk = 1, . . . , n, e la matrice di Householder che annulla gli elementisottodiagonali della colonnak−esima della matrice a cui viene applicata eVk,per k = 1, . . . , n − 2, e la matrice di Householder che annulla gli elementi dipostok+2, . . . , n della rigak−esima della matrice a cui viene applicata.E facile


vedere che, ponendoUB = U1 . . . Un e VB = V1 . . . Vn−2, abbiamo ottenuto lafattorizzazione (2.19). Il costo della trasformazione eO(mn2 − n3/3) se siamointeressati soltanto al calcolo diB; se invece si vogliono calcolare esplicitamenteanche le due matriciVB eUB, il costo aggiuntivo eO(n3/3) eO(mn2 − n3/3)rispettivamente [21].

La bidiagonalizzazione diA e descritta nell’algoritmo 2.3.1, dove si e sceltodi inserire fra i risultati ancheUB eVB. Come nell’analogo algoritmo 1.4.1 per lariduzione di una matrice quadrata in forma di Hessenberg superiore, abbiamo sup-posto di utilizzare matrici di Householder con la scelta (1.9) per i vettoriv che ledefiniscono. In [21] si puo trovare un algoritmo diverso, incui i vettori di House-holder necessari al calcolo delle trasformazioniU1, . . . , Un e V1, . . . , Vn−2, nor-malizzati in modo che la prima componente sia uguale a 1, vengono memorizzatial posto degli elementi diA che risultano azzerati.

Algoritmo 2.3.1. Algoritmo per la riduzione in forma bidiagonale superiore diuna matriceA ∈ R

m×n.Dati: ARisultati: B bidiagonale superiore eUB eVB ortogonali tali cheUT

BAVB = B1.B = A2. UB = Im,VB = In

3. Perk = 1, . . . , n1. v = B(k : m,k)2. v(1) = v(1) + sign(v(1))‖v‖2

3. v = v/‖v‖2

4. B(k : m,k : n) = B(k : m,k : n) − 2v(vTB(k : m,k : n))5. UB(1 : m,k : m) = UB(1 : m,k : m) − 2(UB(1 : m,k : m)v)vT

6. sek ≤ n− 2, allora:v = B(k, k + 1 : n)v(1) = v(1) + sign(v(1))‖v‖2

v = v/‖v‖2

B(k : m,k + 1 : n) = B(k : m,k + 1 : n) − 2(B(k : m,k + 1 : n)v)vT

VB(1 : n, k + 1 : n) = VB(1 : n, k + 1 : n) − 2(VB(1 : n, k + 1 : n)v)vT

Fine sceltaFine ciclo

4. Fine


A =

1.1 2 3.54 −0.5 6−7 0.8 94.1 0.3 1.11 −1 3

.

Applicando adA l’algoritmo precedente in doppia precisione si ottiene la matrice

90 Capitolo 2

B =

13.5285 −1.0485 −0.0000−0.0000 −2.4334 5.9636

0.0000 0.0000 −10.1417−0.0000 0.0000 −0.0000−0.0000 −0.0000 0

,

dove gli elementi “nulli” hanno ordine di grandezza10−15 o inferiore. La normadella matriceUT

BAVB −B risulta5.92e−15. �

2.3.3 Metodo QR

Il metodo QR per calcolare la SVD di una matrice bidiagonale consiste nella ap-plicazione del metodo QR per matrici simmetriche alla matrice tridiagonale sim-metricaT = BTB, senza di fatto formare la matriceT stessa. Vediamo come siapplica il metodo. Data una matriceB ∈ R

n×n bidiagonale superiore

B =

α1 β1 0 0 · · · 00 α2 β2 0 · · · 00 0 α3 β3 · · · 0

. . . . . . . . .. . . αn−1 βn−1

0 0 0 0 0 αn

, (2.20)

si ottiene

T = BTB =

α21 α1β1 0 0 · · · 0

α1β1 α22 + β2

1 α2β2 0 · · · 00 α2β2 α2

3 + β22 α3β3 · · · 0

. . . . . . .. .. . . α2

n−1 + β2n−2 αn−1βn−1

0 0 0 0 αn−1βn−1 α2n + β2

n−1

.

Senza perdita di generalita possiamo supporre che laB sia irriducibile, da cuisegue che anche laT gode di questa proprieta e pertanto, in virtu del corollario1.6.1, ha autovalori tutti semplici.

Un’iterazione del metodo QR per laT consiste nei seguenti passi:

1. Calcolare lo shift, ad esempio l’autovaloreσ di

D =

(

α2n−1 + β2

n−2 αn−1βn−1

αn−1βn−1 α2n + β2

n−1

)

piu vicino aα2n + β2

n−1.


2. Calcolarec = cosθ1 es = sinθ1 tali che(

c s−s c

)T (

α21 − σα1β1

)

=

(

d0

)

,

costruireG1 = G(1, 2, θ1) e calcolareGT1 TG1;

3. Continuare con le successive trasformazioni di GivensG2, . . . , Gn−1 fino aottenere

T = ZTTZ

ancora tridiagonale superiore (Z = G1 · · ·Gn−1).

Una strategia per eseguire l’iterazione senza costruireT fu individuata nel1965 da Golub e Kahan [22]. La matrice bidiagonaleB viene trasformata in unanuova matriceB anch’essa bidiagonale in modo tale cheBT B coincida conT . Siinizia ancora calcolandoG1, che viene poi applicata (a sinistra) aB, dando luogoa una matrice della forma

BG1 =

∗ ∗ 0 0 · · · 0⋄ ∗ ∗ 0 · · · 00 0 ∗ ∗ · · · 0

. .. . . . . . .. . . ∗ ∗

0 0 0 0 0 ∗

.

Si tratta ora di procedere con una strategiabulge-chasingper far scorrere verso ilbasso l’elemento⋄ nato in posizione (2,1). Scegliamo a questo scopo una matricedi GivensU1 tale cheUT

1 BG1 abbia zero al posto (2,1); questo fara nascere unelemento diverso da zero in posizione (1,3) percheUT

1 BG1 ha in generale laforma

UT1 BG1 =

∗ ∗ ⋄ 0 · · · 00 ∗ ∗ 0 · · · 00 0 ∗ ∗ · · · 0

. . . . .. . ... .. ∗ ∗

0 0 0 0 0 ∗

;

per annullare il nuovo⋄ useremo una matrice di GivensV2 = G(2, 3, θ)T conθopportuno, ottenendo

UT1 BG1V2 =

∗ ∗ 0 0 · · · 00 ∗ ∗ 0 · · · 00 ⋄ ∗ ∗ · · · 0

.. . . . . . . .. . . ∗ ∗

0 0 0 0 0 ∗

.

92 Capitolo 2

Possiamo a questo punto procedere con trasformazioni di Givens a destra e asinistra fino ad arrivare a

B = UTn−1 · · ·UT

1 BG1V2 · · ·Vn−1 ≡ UTBV

che e ancora bidiagonale. Inoltre,

BT B = V TBTBV = V TT V ,

dove la prima colonna diV e la stessa diG1 e quindi diZ. Invocando il teoremadel Q implicito, possiamo dire cheBT B coincide con laT che avremmo ottenutoeseguendo un’iterazione del metodo QR sullaT .

L’algoritmo di Golub e Kahan e stato successivamente modificato da Dem-mel e Kahan in [10] nel 1990. Un’attenta analisi dell’erroreha permesso a questidue autori di dimostrare che in presenza di valori singolari“piccoli” e convenienteusare il metodo senza shift perche lo shift puo far perdereaccuratezza. In pocheparole, il motivo di questa maggiore accuratezza e che il metodo senza shift puoessere realizzato senza eseguire sottrazioni, e quindi evitando completamente ilrischio di errori di cancellazione. Demmel e Kahan hanno quindi proposto, stu-diato e sperimentato un algoritmo ibrido in cui si usa il metodo senza shift, salvopassare a quello con shift (che e piu veloce) in base a criteri che facciano riteneredi poterlo fare senza introdurre errori inaccettabili. Il sottoprogramma xBDSQRdi LAPACK, utilizzato dal driver xGESVD, e la funzionesvddi Matlab realizzanol’algoritmo di Demmel e Kahan.

Il sottoprogramma xBDSQR usa il metodo QR quando l’utente richiede ilcalcolo dei valori singolari e dei vettori singolari destrie/o sinistri. Altrimenti, sesi desiderano soltanto i valori singolari, viene utilizzato un metodo diverso, notocon il nome di metododqds. Questo metodo, proposto da Fernando e Parlett in[18] nel 1994 e riassunto in [8], calcola a partire daB0 = B una successione{Bi}di matrici bidiagonali, tutte con gli stessi valori singolari; la successione convergead una matrice diagonale che ha quindi sulla diagonale i valori singolari diB.Fondamentalmente, si dimostra che eseguire due iterazionidi questo metodo eequivalente a eseguire un’iterazione del metodo QR con shift (per un’opportunascelta dello shift); da questo segue che la convergenza e molto veloce. Non solo,ma il metododqdsgarantisce alta accuratezza per tutti i valori singolari, indipen-dentemente dalla loro grandezza, perche puo venire realizzato senza necessita dieseguire sottrazioni. Il fatto che questo metodo non puo essere utilizzato per cal-colare i vettori singolari dipende dal fatto che il passaggio daBi a Bi+1 nonutilizza trasformazioni ortogonali.

2.3.4 Metodo di Jacobi

Nel paragrafo 1.6.5 abbiamo visto il metodo di Jacobi per calcolare autovalorie autovettori di matrici simmetriche. Le caratteristiche fondamentali di questometodo sono due, una positiva e l’altra negativa. La possibilita di calcolare con


alta accuratezza gli autovalori a prescindere dalla loro grandezza rappresenta illato positivo, ed e dovuta al fatto che la matrice non viene inizialmente tridiago-nalizzata. L’aspetto negativo del metodo e la sua bassa velocita di convergenza,soprattutto se confrontata con quella dei metodi basati sulla tridiagonalizzazione.

Nel campo dei metodi per calcolare la SVD si presenta esattamente la stes-sa situazione. Il metodo QR e le sue varianti sono molto efficienti perche basatisulla bidiagonalizzazione diA, ma la trasformazione iniziale daA aB puo esserepenalizzante dal punto di vista dell’accuratezza, rendendo con cio questi metodinon adatti a tutte quelle applicazioni in cui i valori singolari piu piccoli sono par-ticolarmente importanti. Per far capire meglio il problema, riportiamo un esempio“estremo” tratto da [8].

Esempio 2.3.2.Siaη una costante positiva. Consideriamo la matrice

A =

η 1 1 1η η 0 0η 0 η 0η 0 0 η

,

i cui valori singolari sono√

3,√

3η, η, η. Supponiamo di bidiagonalizzareAtramite trasformazioni di Householder, usando ad esempio l’algoritmo 2.3.1. Sepotessimo operare in precisione infinita, al primo passaggio di questa trasfor-mazione si otterrebbe la matrice

U1A =

−2η −0.5 − η2 −0.5 − η

2 −0.5 − η2

0 −0.5 + 5η6 −0.5 − η

6 −0.5 − η6

0 −0.5 − η6 −0.5 + 5η

6 −0.5 − η6

0 −0.5 − η6 −0.5 − η

6 −0.5 + 5η6

.

D’altra parte, dobbiamo lavorare in aritmetica floating-point. Supponiamo che laprecisione di macchina siaǫm ≃ 10−16 e cheη sia molto piccolo, per esempioη = 10−20; allora si ottiene

fl(U1A) =

−2η −0.5 −0.5 −0.50 −0.5 −0.5 −0.50 −0.5 −0.5 −0.50 −0.5 −0.5 −0.5

e questa matrice ha rango 2 e valori singolari√

3,√

3η, 0, 0. Quindi, gia conquesto primo passaggio ci siamo giocati la possibilita di calcolare accuratamentegli ultimi due valori singolari diA. Nella realta, proseguendo con la bidiagonaliz-zazione in precisione finita, arriveremo a una matriceB i cui ultimi due valori sin-golari non saranno esattamente uguali a zero, ma saranno “piccoli” semplicementeper effetto degli errori di arrotondamento. A titolo di esempio, con il comandosvd(A) di Matlab che utilizza il metodo QR passando dalla bidiagonalizzazione,si ottengono i seguenti valori (riportiamo soltanto le prime cifre per comodita):

94 Capitolo 2

1.7321e− 00, 6.7987e− 17, 6.2804e− 36, 0. Il primo di questi valori e accuratoa precisione di macchina, mentre gli altri tre sono completamente sbagliati.

Osserviamo che le stesse considerazioni potremmo farle seA avesse la forma

A =

η α α αη η 0 0η 0 η 0η 0 0 η

conα eη costanti positive, i cui valori singolari sono quelli di prima moltiplicatiperα. Se fosseη << 1 + αǫm incorreremmo nei soliti problemi di prima. �

Passiamo ora a descrivere il metodo di Jacobi, supponendo per semplicita divoler calcolare la SVD ridotta. Ripercorriamo il metodo di Jacobi per la diago-nalizzazione diATA. In linea di principio si dovrebbe costruire la successione

{

Λ0 = ATAΛk = GT

k Λk−1Gk,

doveGk e una matrice di GivensG(p, q, θ) scelta in modo che la matriceΛk

abbia gli elementi di posto(p, q) e (q, p) uguali a zero. La successione{Λk}convergera alla matrice diagonaleΛ degli autovalori diATA, in cui si potranno aposteriori ordinare gli autovalori in senso decrescente, ela successione{Vk}, conVk = G1G2 · · ·Gk, tendera alla matrice dei vettori singolari destriV .

Guardiamo da vicino la successione. Perk = 1, abbiamoΛ1 = GT1 A

TAG1;se poniamoA1 = AG1, possiamo scrivereΛ1 = AT

1A1. Alla successiva itera-zione, abbiamoΛ2 = GT

2AT1A1G2; se poniamoA2 = A1G2, possiamo scrivere

Λ2 = AT2 A2. E cosı via. Queste notazioni suggeriscono allora di non creare la

successione{Λk}, ma la{Ak}, ossia

{

A0 = AAk = Ak−1Gk,

(2.21)

che pertanto tendera aΣ. In virtu della definizione (2.21), in letteratura il metodoviene chiamato metodo di Jacobione-sidede fu proposto gia nel 1958 da M.R.Hestenes in [25].

Un paio di osservazioni di carattere pratico sono indispensabili, relative aquei punti dell’algoritmo nei quali potrebbe nascere l’esigenza di calcolare la ma-triceAT

k−1Ak−1, cosa che ovviamente non possiamo fare. Come prima cosa, percreare la matriceAk daAk−1 dobbiamo scegliereGk = G(p, q, θ) in modo daannullare l’elemento di posto(p, q) di AT

k−1Ak−1. Ricordando la formula (1.50),vediamo che a questo scopo non occorre calcolare tutta la matrice ma soltanto isuoi elementi di posto(p, q), (p, p) e (q, q). La seconda osservazione riguardala variante del metodo di Jacobi da utilizzare. La variante classica non va beneperche imporrebbe il calcolo di tutta la matriceAT

k−1Ak−1 ad ogni iterazione per


individuare gli indici(p, q) dell’elemento extradiagonale piu grande in valore as-soluto; dobbiamo quindi necessariamente usare la versioneciclica del metodo. Aquesto proposito, ricordiamo che in un algoritmo di Jacobi ciclico si esegue larotazioneG(p, q, θ) solo se l’elemento di indici(p, q) non e troppo piccolo: nelcontesto della SVD, Demmel in [8] propone il criterio

|bpq| > ǫm√

bppbqq.

Demmel e altri colleghi si sono occupati negli anni ’90 di individuare clas-si di matrici per i quali il metodo di Jacobi one-sided garantisce alta accuratezzanel calcolo di tutti i valori singolari. Una classe, ad esempio, e quella delle ma-trici A quadrate che possono essere fattorizzate nella formaA = DX, conDdiagonale non singolare eX ben condizionata. A questa classe appartiene anchela matrice dell’esempio 2.3.2, che infatti puo essere ottenuta come prodotto diD = diag(1, η, η, η) e

X =

η 1 1 11 1 0 01 0 1 01 0 0 1

.

Perη = 1e−20, il numero di condizionamento diX, calcolato con la funzioneconddi Matlab, e 2.3028.

Il metodo di Jacobi one-sided e molto accurato, spesso piudei metodi basatisul QR, e contemporaneamente eredita la lentezza del metododi Jacobi per ilproblema degli autovalori. Ma gli analisti numerici sono spesso ostinati, e alcuniricercatori hanno recentemente indirizzato i loro sforzi aimmaginare varianti delmetodo che infrangessero la barriera della lentezza. Gli sforzi sono stati premiatie nel 2008 e stato pubblicato un articolo da Z. Drmac e K. Veselic [14] in cui sipresenta una variante del metodo di Jacobi che eredita tuttele buone proprieta diaccuratezza del metodo originale e puo superare il metodo QR dal punto di vistadell’efficienza. Il risultato e stato possibile grazie all’utilizzo di precondizionatori,uno strumento potente in molti campi dell’analisi numerica. Ma questa e un’altrastoria....

96 Capitolo 2

AVVERTENZANonostante tutta l’attenzione che ho cercato di mettere nella stesura di questi

appunti e le riletture mie e di volonterosi colleghi, sono consapevole che si trattaancora sostanzialmente di una prima stesura. Questo implica che alcuni errori,forse non solo di stampa, sono sicuramente rimasti e che alcune parti sono ancorascritte in modo non completamente soddisfacente. Ovviamente me ne assumotutta la responsabilita e sono grata fin da ora a chiunque mi aiutera a migliorare iltutto segnalandomi qualunque osservazione.Il mio indirizzo e-mail e:[email protected].

Maria Grazia Gasparo

Firenze, Dicembre 2010

Bibliografia

[1] E. Anderson e altri,LAPACK User’s guide, third edition, SIAM, 1999[2] L. Autonne, Sur les matrices hypohermitiennes et les unitairs, Comptes

rendus de l’Academie Sciences 156 (1913), 858-860[3] J.L. Barlow,Error analysis of update methods for the symmetric eigenvalue

problem, SIAM J. MAtrix Anal. Appl. 14 (1993), 598-618[4] W. Barth, R.S. Martin, J.H. Wilkinson,Calculation of the eigenvalues of

a symmetric tridiagonal matrix by the method of bisection, Numer. Math.9 (1967), 386-393, ripubblicato in J.H. Wilkinson e C. Reinsch, Linearalgebra, Handbook for automatic computation, Vol.II, Springer-Verlag, 1971

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

[6] S. Chandrasekaran, I.C.F. Ipsen,On rank-revealing factorizations, SIAM J.Matrix Anal. Appl., 15 (1994), 592-622

[7] J.J.M. Cuppen,A Divide and Conquer method for the symmetric tridiagonaleigenproblem, Numer. Math. 36 (1981),177-195

[8] J.W. Demmel,Applied numerical linear algebra, SIAM, 1997[9] J.W. Demmel, I. Dhillon, H. Ren,On the correctness of some bisection-like

parallel eigenvalue algorithms in floating-point arithmetic, Electronic Trans.Numer. Anal. (ETNA), 3 (1995) 116-149

[10] J. W. Demmel, W. Kahan,Accurate singular values of bidiagonal matrices,SAM J. Sci. Stat. Comput. 11 (1990) 873-912

[11] J.W. Demmel, O.A. Marques, B.N. Parlett, C. Vomel, Performance andaccuracy of LAPACK’s symmetric tridiagonal eigensolvers, SIAM J. Sci.Comput. 30 (2008), 1508-1526

[12] J.W. Demmel, K. Veselic,Jacobi’s method is more accurate than QR, SIAMJ. Matrix Anal. Appl. 13 (1992), 1204-1245

[13] J.J. Dongarra, D.C. Sorensen,A fully parallel algorithm for the symmetriceigenproblem, SIAM J. Sci. Statist. Comput. 8 (1987), 139-154

[14] Z. Drmac, K. Veselic,New fast and accurate Jacobi SVD algorithm. I, SIAMJ. Matrix Anal. Appl. 29 (2008), 1322-1342

[15] C. Eckart, G. Young,The approximation of one matrix by another of lowerrank, Psychometrika 1 (1936), 211-218

[16] C. Eckart, G. Young,A principal axis transformation for non-hermitianmatrices, Bull. Amer. Math. Soc. 45 (1939), 118-121

[17] A.N. Elmachtaub, C.F. Van Loan,From random polygon to ellipse, SIAMReview 52 (2010), 151-170

98 Bibliografia

[18] V. Fernando, B. Parlett,Accurate singular values and differential qdalgorithms, Numer. Math. 67 (1994) 191-229

[19] F. Fontanella, A. Pasquali,Calcolo numerico. Metodi ed algoritmi, Pitagoraeditrice, 1982

[20] J.W. Givens,Numerical computation of the characteristic values of a realsymmetric matrix, Oak Ridge National Laboratory, ORNL-1574 (1954)

[21] G.H. Golub, C.F. Van Loan,Matrix computations, terza ed., Johns HopkinsUniversity Press, Baltimore, 1996

[22] G.H. Golub, W. Kahan,Calculating the singular values and pseudoinversesof a matrix, SIAM J. Numer. Anal. Ser.B, 2 (1965) 205-224

[23] G.H. Golub, F. Uhlig,The QR–algorithm: 50 years later its genesis byJohn Francis and Vera Kublanovskaya, and subsequent developments, IMAJournal of Numerical Analysis, 29 (2009), 467-485

[24] P. Henrici,On the speed of convergence of cyclic and quasicyclic Jacobimethod for computing eigenvalues of hermitian matrices, J. Soc. Ind. Appl.Math. 6 (1958), 144-162

[25] M.R. Hestenes,Inversion of matrices by biorthogonalization and relatedresults, J. Soc. Indust. Appl. Math. 6 (1958), 51-90

[26] R.A. Horn, C.R. Johnson,Matrix analysis, Cambridge University Press,1985

[27] D.A. Huckaby, T.F. Chan,Stewart’s pivoted QLP decomposition for low-rank matrices, Numer. Lin. Algebra with Appl. 12 (2005), 153-159

[28] C.G.J. Jacobi,Uber ein leichtes Verfahren, die in der Theorie derSacularstorungen vorkommenden Gleichungen numerische aufzulosen, J.reine u. angew. Math., 30 (1846), 51-95

[29] R.S. Martin, G. Peters, J.H. Wilkinson,The QR algorithm for real Hessen-berg matrices, in J.H. Wilkinson e C. Reinsch,Linear algebra, Handbookfor automatic computation, Vol.II, Springer-Verlag, 1971, pag.359-371

[30] C.L. Muntz, Solution directe de l’equation seculaire et de quelquesproblemes analogues trascendants, Comptes Rendus de l’Academie desSciences de Paris, 156 (1913), 43-46

[31] A.M. Ostrowski, On the convergence of the Rayleigh quotient iterationfor the computation of the characteristic roots and vectors.I, Archives forRational Mechanics and Analysis, 1 (1957), 233-241

[32] B.N. Parlett,The symmetric eigenvalue problem, Prentice-Hall, EnglewoodCliffs, NJ, 1980, ristampato anche come Classics in AppliedMathematics20, SIAM, Philadelphia, 1997

[33] G. Peters, J.H. Wilkinson,Inverse iteration, ill–conditioned equations andNewton’s method, SIAM Review, 3 (1979), 339-360

[34] J. Rutter,A serial implementation of Cuppen’s Divide and conquer algo-rithm for the symmetric tridiagonal eigenproblem, Computer Science Di-vision Report UCB/CSD 94/799, University of California, Berkeley, 1994(anche LAPACK Working Note 69)

[35] Y. Saad,Numerical methods for large eigenvalue problems, ManchesterUniversity press, 1992

Bibliografia 99

[36] G.W. Stewart,On the early history of the Singular Value Decomposition,SIAM Review 35 (1993), 551-566

[37] G.W. StewartMatrix Algorithms Vol. I: Basic Decompositions, SIAM, 1998[38] G.W. StewartMatrix Algorithms Vol. II: Eigensystems, SIAM, 1998[39] J. Stoer, R. Bulirsch,Introduzione all’analisi numerica, Vol. 2, trad. di G.C.

Barozzi, Zanichelli ed., 1975[40] L.N. Trefethen, D. Bau III,Numerical Linear Algebra, SIAM, 1997[41] H.P.M. van Kempen,On the quadratic convergence of the special cyclic

Jacobi method, Numer. Math. 9 (1966),19-22[42] D.S. Watkins,The QR algorithm revisited, SIAM Review 50 (2008), 133-

145[43] H. Wielandt, Beitrage zur mathematischen Behandlung komplexer Eigen-

wertprobleme, Teil V: Bestimmung Hoherer Eigenwerte durch gebrocheneIteration, Bericht B 44/J/37, Aerodynamische Versuchsanstalt Gottingen,Germany, 1944

[44] J.H. Wilkinson, Note on the quadratic convergence of the cyclic Jacobiprocess, Numer. Math. 4 (1962), 296-300

[45] J.H. Wilkinson,The algebraic eigenvalue problem, Clarendon Press, Oxford,1965

[46] J.H. Wilkinson,Global convergence of tridiagonal QR algorithm eith originshifts, Lin. Alg. and its Applic.,1 (1968),409-420

Maria Grazia Gasparo Metodi numerici per il calcolo di ... · L’algoritmo usato da Google e`...

Documents

Transcript of Maria Grazia Gasparo Metodi numerici per il calcolo di ... · L’algoritmo usato da Google e`...