ARGOMENTI DEL CORSO CALCOLO NUMERICOvvw.web.cs.unibo.it/wiki/images/3/3a/CC-2009_2010-05.pdf1....

ARGOMENTI DEL CORSO

CALCOLO NUMERICO

A.A. 2009/10

A = L · U

L=

1

ℓ2,1 1...

. . .

ℓn,1 . . . ℓn,n−1 1

U =

u1,1 . . . u1,n

u2,2. . .

...

un,n

Algebra Lineare Numerica

Giulio Casciola

(novembre 2004, rivista e corretta novembre 2009)

Indice

1 Sistemi Lineari: metodi diretti 11.1 Fattorizzazione LU . . . . . . . . . . . . . . . . . . . . . . 3

1.1.1 Sostituzione in Avanti . . . . . . . . . . . . . . . . 31.1.2 Sostituzione all’indietro . . . . . . . . . . . . . . . . 41.1.3 Metodo di Gauss . . . . . . . . . . . . . . . . . . . 61.1.4 Calcolo di A−1 . . . . . . . . . . . . . . . . . . . . 91.1.5 Calcolo del det(A) . . . . . . . . . . . . . . . . . . 101.1.6 Metodo di Gauss con scambio delle righe . . . . . . 101.1.7 Esistenza della Fattorizzazione P A = L U . . . . . 121.1.8 Stabilita della Fattorizzazione L U . . . . . . . . . 131.1.9 Metodo di Gauss con scambio delle righe e perno

massimo . . . . . . . . . . . . . . . . . . . . . . . . 141.2 Condizionamento del Problema Ax = b . . . . . . . . . . . 14

1.2.1 Richiami sul concetto di norma . . . . . . . . . . . 151.2.2 Errore Inerente . . . . . . . . . . . . . . . . . . . . 17

1.3 Fattorizzazione QR . . . . . . . . . . . . . . . . . . . . . . 191.3.1 Matrici elementari di Householder . . . . . . . . . . 201.3.2 Metodo di Householder . . . . . . . . . . . . . . . . 221.3.3 Implementazione del metodo di Householder . . . . 251.3.4 Costo Computazionale per risolvere Ax = b tramite

Q R . . . . . . . . . . . . . . . . . . . . . . . . . . 251.3.5 Stabilita Numerica della Fattorizzazione Q R . . . 26

2 Autovalori e Autovettori 272.1 Proprieta degli Autovalori . . . . . . . . . . . . . . . . . . 292.2 Proprieta degli Autovettori . . . . . . . . . . . . . . . . . . 302.3 Similitudine fra Matrici . . . . . . . . . . . . . . . . . . . . 31

3 Sistemi lineari: metodi iterativi 353.1 Decomposizione della matrice . . . . . . . . . . . . . . . . 363.2 Controllo della Convergenza . . . . . . . . . . . . . . . . . 39

3

4 INDICE

3.3 Test di Arresto . . . . . . . . . . . . . . . . . . . . . . . . 403.4 Metodi di Jacobi e Gauss-Seidel . . . . . . . . . . . . . . . 41

4 Metodi per Autovalori e Autovettori 454.1 Riduzione di una matrice in forma di Hessenberg . . . . . 474.2 Metodo Q R per il calcolo degli autovalori . . . . . . . . . 47

4.2.1 Algoritmo di base . . . . . . . . . . . . . . . . . . . 484.2.2 Risultati di convergenza . . . . . . . . . . . . . . . 484.2.3 Costo computazionale e stabilita . . . . . . . . . . 50

4.3 Condizionamento del calcolo degli autovalori . . . . . . . . 50

5 Il problema dei Minimi Quadrati 535.1 Le equazioni normali . . . . . . . . . . . . . . . . . . . . . 535.2 Metodo QR per i minimi quadrati . . . . . . . . . . . . . . 54

Bibliografia 57

Capitolo 1

Sistemi Lineari: metodi diretti

Uno dei problemi piu frequenti nel calcolo scientifico e la soluzione di unsistema lineare; se questo e composto da tante equazioni quante incognitesi dice sistema quadrato o normale. In forma matriciale puo essere scrittocome

Ax = b

dove A e una data matrice di ordine n × n, b e un dato vettore colonnacon n elementi ed x e il vettore delle incognite.In algebra lineare si studiano metodi per risolvere sistemi lineari non sin-golari. Un metodo noto e quello di Cramer (od anche regola di Cramer)nel quale ogni componente della soluzione e espressa come quoziente di undeterminante sempre diverso a numeratore e di uno stesso determinante adenominatore, cosı che la soluzione del sistema lineare si riduce al calcolodi n + 1 determinanti di ordine n. Se si cerca di risolvere un sistema di20 equazioni con la regola di Cramer, sarebbe necessario calcolare 21 de-terminanti di ordine 20. Come e noto per calcolare il determinante di unamatrice n× n, per esempio con la formula di Laplace o quella di Leibniz,servono (n−1)n! moltiplicazioni; nel nostro esempio numerico sara 19 ·20!e quindi l’intero sistema lineare comportera 19 · 20! · 21 moltiplicazioni,piu un ugual numero di addizioni. Su un normale PC, oggi si possono farecirca 109 moltiplicazioni al secondo (Giga FLOPS FLoating point OPera-tion per Second), cosı che solo le moltiplicazioni richiederanno circa 30782anni, sempre che non manchi la corrente durante il calcolo.In algebra lineare si insegna anche che la soluzione di Ax = b puo esserescritta come x = A−1b, dove A−1 e l’inversa di A. Comunque nella mag-gioranza dei problemi pratici, cio non e necessario ed eventualmente risul-ta inopportuno calcolare A−1. Con un esempio estremo, ma illustrativo,

1

2 CAPITOLO 1. SISTEMI LINEARI: METODI DIRETTI

consideriamo un sistema di appena una equazione, come

7x = 21.

Il miglior modo per risolvere tale sistema e con la divisione

x =21

7= 3.

L’uso della matrice inversa porterebbe invece a

x = 7−1 · 21 = 0.142857 . . . · 21 = 2.99997 . . . .

L’inversa richiede piu operazioni, una divisione ed una moltiplicazioneinvece di appena una divisione. E il maggior numero di operazioni checi induce ad evitare il calcolo dell’inversa. Ci concentreremo quindi sullasoluzione diretta di sistemi piuttosto che sul calcolo dell’inversa.

Vedremo essenzialmente due metodi diretti per risolvere un sistemalineare A x = b. Entrambi i metodi si basano su una fattorizzazione dellamatrice A, ossia costruiscono due matrici il cui prodotto sia uguale adA:

1. Fattorizzazione LU (L (Low) matrice triangolare inferiore e diago-nale unitaria, U (Up) matrice triangolare superiore)

A = L U

2. Fattorizzazione QR (Q matrice ortogonale, R (Right) matrice trian-golare superiore)

A = Q R

Lo scopo e quello di ricondurre la soluzione del sistema (pieno e/o senzastruttura particolare) Ax = b alla soluzione di uno o due sistemi (sparsie/o con particolare struttura); la soluzione del sistema si ottiene in duefasi:

1. fattorizzazione di A;

2. soluzione di sistemi lineari di forma triangolare, le cui caratteristichedi soluzione sono piu semplici e meno costose (rispetto a quelle delsistema originario).

1.1. FATTORIZZAZIONE LU 3

1.1 Fattorizzazione LU

In questa sezione si studia la possibilita di fattorizzare la matrice A nel pro-dotto di una matrice L triangolare inferiore per una matrice U triangolaresuperiore:

An×n = Ln×n · Un×n

con

L =

1ℓ2,1 1ℓ3,1 1

.... . .

ℓn,1 . . . ℓn,n−1 1

U =

u1,1 . . . u1,n

u2,2

u3,3

. . ....

un,n

.

L’esistenza di una tale fattorizzazione renderebbe facile la determinazionedella soluzione del sistema normale non singolare, infatti si avrebbe:

LUx = b

e ponendo Ux = y si potrebbe risolvere il sistema

Ly = b

per sostituzione in avanti (forward-substitution); determinato y si risolveil sistema

Ux = y

per sostituzione all’indietro (backward-substitution).

1.1.1 Sostituzione in Avanti

Il primo sistema da risolvere e

L y = b

in cui L e una matrice triangolare inferiore.

• La soluzione si ottiene con la sostituzione in avanti, mediante l’algo-ritmo che segue

y(1) = b(1)/L(1, 1)

per i = 2, ..., n

y(i) = b(i)− L(i, 1 : i− 1) ∗ y(1 : i− 1)


y(i) = y(i)/L(i, i)

La struttura dell’algoritmo e triangolare. Per generare y(i) si com-piono 2i − 1 operazioni e per la precisione i − 1 moltiplicazioni, 1divisione e i− 1 sottrazioni.

• Il costo computazionale della sostituzione in avanti e dato da

1 +n∑

i=2

(2i− 1) =n∑

i=1

(2i− 1) = 2n∑

i=1

i− n = n(n + 1)− n

• Se L e una matrice triangolare inferiore con elementi diagonali uni-tari, allora nell’algoritmo non ci sono divisioni:

y(1) = b(1)

per i = 2, ..., n

y(i) = b(i)− L(i, 1 : i− 1) ∗ y(1 : i− 1)

ed il costo computazionale diventa

n∑

i=2

2(i− 1) =n∑

i=1

2(i− 1) = n(n + 1)− 2n

Esempio 1.1 Sostituzione in avanti.

ℓ1,1 0 0ℓ2,1 ℓ2,2 0ℓ3,1 ℓ3,2 ℓ3,3

y1

y2

y3

=

b1

b2

b3

y1 = b1/ℓ1,1 1 operazione↓y2 = (b2 − ℓ2,1 ∗ y1)/ℓ2,2 3 operazioni↓y3 = (b3 − ℓ3,1 ∗ y1 − ℓ3,2 ∗ y2)/ℓ3,3 5 operazioni

Costo computazionale: n2 = 32 = 9 operazioni.

1.1.2 Sostituzione all’indietro

Il secondo sistema da risolvere e

U x = y

in cui U e una matrice triangolare superiore.


• La soluzione si ottiene con la sostituzione all’indietro, mediantel’algoritmo che segue

x(n) = y(n)/U(n, n)

per i = n− 1, ..., 1

x(i) = y(i)− U(i, i + 1 : n) ∗ x(i + 1 : n)

x(i) = x(i)/U(i, i)

La struttura dell’algoritmo e triangolare. Per generare x(i) si com-piono 2i − 1 operazioni e per la precisione i − 1 moltiplicazioni, 1divisione e i− 1 sottrazioni.

• Il costo computazionale della sostituzione all’indietro e dato da

1 +n∑

i=2

(2i− 1) = ... = n2

Esempio 1.2 Sostituzione all’indietro

u1,1 u1,2 u1,3

0 u2,2 u2,3

0 0 u3,3

x1

x2

x3

=

y1

y2

y3

x3 = y3/u3,3 1 operazione↓x2 = (y2 − u2,3 ∗ x3)/u2,2 3 operazioni↓x1 = (y1 − u1,2 ∗ x2 − u1,3 ∗ x3)/u1,1 5 operazioni

Costo computazionale: n2 = 32 = 9 operazioni.

Riassumendo: determinata una fattorizzazione LU di A, il vettore solu-zione x si puo ottenere facilmente risolvendo due sistemi triangolari, unoinferiore e l’altro superiore.

Ci si concentrera allora sulla fattorizzazione LU di A; si osserva chela sola ipotesi che A non sia singolare non garantisce l’esistenza di unafattorizzazione LU , infatti vale il seguente teorema.

Teorema 1.1 Se i minori principali di ordine k di A per k = 1, . . . , n−1sono diversi da zero, allora esiste una ed una sola fattorizzazione LU diA.


1.1.3 Metodo di Gauss

Procediamo in un caso semplice ad illustrare come funziona il metodo diGauss per fattorizzare LU una matrice An×n (il procedimento in oggettopuo essere applicato anche a matrici rettangolari).

Sia

A =

a1,1 a1,2 a1,3

a2,1 a2,2 a2,3

a3,1 a3,2 a3,3

con

a1,1 6= 0 e det

(

a1,1 a1,2

a2,1 a2,2

)

= a1,1a2,2 − a1,2a2,1 6= 0.

Consideriamo la matrice L1 siffatta:

L1 =

1 0 0

−a2,1a1,1

1 0

−a3,1a1,1

0 1

(si noti la necessita che a1,1 6= 0; a1,1 e detto perno o pivot) allora sara:

L1A =

a1,1 a1,2 a1,3

0 a(1)2,2 a

(1)2,3

0 a(1)3,2 a

(1)3,3

con a(1)i,j = ai,j − ai,1

a1,1a1,j, i, j = 2, 3.

Consideriamo ora la matrice L2 siffatta:

L2 =

1 0 00 1 0

0 −a(1)3,2

a(1)2,2

1

(si noti la necessita che a(1)2,2 =

a1,1a2,2 − a1,2a2,1a1,1

6= 0, ossia che a1,1a2,2 −a1,2a2,1 6= 0; a

(1)2,2 e detto perno o pivot) allora sara:

L2(L1A) =

a1,1 a1,2 a1,3

0 a(1)2,2 a

(1)2,3

0 0 a(2)3,3


con a(2)3,3 = a

(1)3,3 −

a(1)3,2

a(1)2,2

a(1)2,3.

Chiamiamo U la L2L1A e notiamo che la U e del tipo cercato. Si noti chela L1 e la L2 sono non singolari, allora

A = L−11 L−1

2 U.

Chi sono L−11 e L−1

2 e chi e L−11 L−1

2 ? Saranno:

L−11 =

1 0 0a2,1a1,1

1 0a3,1a1,1

0 1

L−12 =

1 0 00 1 0

0a

(1)3,2

a(1)2,2

1

L−11 L−1

2 =

1 0 0a2,1a1,1

1 0

a3,1a1,1

a(1)3,2

a(1)2,2

1

e se chiamiamo quest’ultima L avremo che A = L U come volevamo.U ←− A; L←− In;

per k = 1, .., n− 1

per j = k + 1, .., n

ℓ(j, k) = u(j, k)/u(k, k)

u(j, k + 1 : n) = u(j, k + 1 : n)− ℓ(j, k) ∗ u(k, k + 1 : n)

Per determinare la fattorizzazione suddetta sono necessarie n(n + 1)/2moltiplicazioni/divisioni floating point per il calcolo degli elementi di L e(n−1)2+(n−2)2+. . .+22+1 = (n−1)n(2n−1)/6 moltiplicazioni/divisionifloating point per gli elemnti di U . La complessita computazionale, intermini di operazioni di moltiplicazione e divisione, e quindi:

(n− 1)n(n + 1)

3=

n3 − n2

3≈ 1

3n3.

Esempio 1.3 Fattorizzazione A = L U di Gauss.

Sia Ax = b con A =

2 1 04 5 26 15 12

e b =

212

.


L1 deve essere triangolare inferiore e tale da rendere nulli gli elementidella prima colonna di A sotto l’elemento a11 = 2.

L1 =

1 0 0−2 1 0−3 0 1

L1 A =

1 0 0−2 1 0−3 0 1

·

2 1 04 5 26 15 12

=

2 1 00 3 20 12 12

L2 deve essere triangolare inferiore e tale da rendere nulli gli elementidella seconda colonna di L1 A sotto ad a

(1)2,2 = 3.

L2 =

1 0 00 1 00 −4 1

L2 (L1 A) =

1 0 00 1 00 −4 1

·

2 1 00 3 20 12 12

=

2 1 00 3 20 0 4

= U

Allora

L2 L1 A = U,

A = L−11 L−1

2 U = L U

e risulta

L = L−11 L−1

2 =

1 0 02 1 03 4 1

.

Fattorizzata la matrice, si procede alla soluzione dei sistemi triangolari

Ly = b e Ux = y.

1 0 02 1 03 4 1

y1

y2

y3

=

212

y1 = 2;

(

1 04 1

) (

y2

y3

)

=

(

−3−4

)

y2 = −3;

(1)(y3) = (8) y3 = 8;


2 1 00 3 20 0 4

x1

x2

x3

=

2−38

x3 = 2;

(

2 10 4

) (

x1

x2

)

=

(

2−7

)

x2 = −7/3;

(2)(x1) = (13/3) x1 = 13/6.

1.1.4 Calcolo di A−1

Sia A ∈ IRn×n, invertibile; vogliamo calcolare A−1, risolvendo il sistema

A X = In, X = [x1 | x2 | ...... | xn], In = [e1 | e2 | ...... | en]

In e la matrice identita di dimensione n; i vettori e1, ... en sono la basecanonica di IRn. X ∈ IRn×n e la matrice incognita.

Risolvere A X = In, per ottenere X = A−1 In, e equivalente a risolvere nsistemi lineari, tutti con la stessa matrice A dei coefficienti

A x1 = e1

A x2 = e2...

A xn = en

Si fattorizzi A = L U una sola volta, e poi si risolvano 2 n sistemi linearidi forma triangolare

(i) A = L U 13n3 − 1

3n2

(ii)

Ly1 = e1 n2 − n

Ux1 = y1 n2

...

...

Lyn = en n2 − n

Uxn = yn n2


1.1.5 Calcolo del det(A)

Fattorizzare A = L U ci permette di calcolare con costo computazionalebasso il determinante della matrice A; vale

det(A) = det(L) det(U) = det(U) =n∏

i=1

ui,i

Abbiamo utilizzato il Teorema di Binet, che det(L) = 1 e che il determi-nante di una matrice triangolare e dato dal prodotto degli elementi delladiagonale.

1.1.6 Metodo di Gauss con scambio delle righe

Sia dato il sistema lineare(

0 31 2

) (

x1

x2

)

=

(

50

)

;

essendo a1,1 = 0, non si puo applicare l’algoritmo di fattorizzazione diGauss alla matrice dei coefficienti, ma si osserva che questo sistema eequivalente a

(

1 20 3

) (

x1

x2

)

=

(

05

)

;

ed ora a1,1 6= 0.Allora se A e non singolare e sempre possibile trovare una matrice dipermutazione P per cui P A sia fattorizzabile L U , cioe

P A = L U.

Matrice di Permutazione: si ottiene dalla matrice identita In permu-tando le righe. Applichiamo P ad una matrice A ∈ IRn×n :

P A −→ permuta righe di A

A P −→ permuta colonne di A

Ad esempio:

A =

1 2 34 5 67 8 9

P =

1 0 00 0 10 1 0

P A =

1 2 37 8 94 5 6

A P =

1 3 24 6 57 9 8

.


Ovviamente non e pensabile determinare in qualche modo una matrice Popportuna, premoltiplicarla per A (scambio di righe e quindi ordine delleequazioni) e poi applicare l’algoritmo visto, bensı si dovra determinare Pcostruttivamente. Questo e cio che fa l’algoritmo di eliminazione di Gausscon scambio delle righe (pivoting parziale). Si illusta il metodo con unesempio. Sia data la matrice

A = A0 =

4 −8 22 −4 61 −1 3

;

il pivot a(0)i,1 e 4 e quindi P1 = I e

L1 =

1 0 0−1/2 1 0−1/4 0 1

;

si costruisce

A1 = L1 P1 A0 =

4 −8 20 0 50 1 5/2

.

Il pivot a(1)i,2 e 1; le matrici P2 ed L2 saranno:

P2 =

1 0 00 0 10 1 0

L2 =

1 0 00 1 00 0 1

.

Si costruisce

A2 = L2 P2 A1 =

4 −8 20 1 5/20 0 5

= U.

In definitiva eL2 P2 L1 P1 A = U.

Si noti cheL2 P2 L1 P1 = L2 P2 L1 P−1

2 P2P1

e se poniamo T = L2 P2 L1 P−12 si ha che T e triangolare inferiore e quindi

T P2 P1 A = U

e ponendoL−1 = T e P = P2 P1


si ha

P A = L U.

Completando l’esempio si ha

L = P2 L−11 P−1

2 L−12 =

1 0 01/4 1 01/2 0 1

ed infatti

L U =

1 0 01/4 1 01/2 0 1

·

4 −8 20 1 5/20 0 5

=

4 −8 21 −1 32 −4 6

= P A.

Nell’implementare tale metodo ci si accorge immediatamente che non epossibile perdere l’informazione P , infatti se fattorizziamo P A il sistemada risolvere sara

P A x = P b.

1.1.7 Esistenza della Fattorizzazione P A = L U

Che cosa possiamo fare se dobbiamo considerare una matrice A per laquale il metodo di fattorizzazione semplice di Gauss fallisce? Arrichiamoil metodo con un procedimento di pivotaggio. Vale infatti il seguenteteorema.

Teorema 1.2 (Esistenza di P A = L U) Per ogni matrice A ∈ IRn×n

non singolare, esiste una matrice di permutazione P tale che

P A = L U

con L triangolare inferiore con elementi unitari e U triangolare superiore.

Il teorema enunciato assicura che, se la fattorizzazione di Gauss diuna matrice A non singolare sembra fallire (perche ad un certo passo k

si ottiene un perno a(k)k,k ≈ 0), e possibile superare questo passo andando

a cercare (nella colonna k−esima della matrice Ak corrente) un nuovo

perno a(k)i,k 6= 0 ed introducendo una permutazione di righe in modo che


a(k)i,k occupi la posizione di a

(k)kk

Lk ... L1 A =

. . ....

. . . a(k)k,k = 0 . . . ← riga k

a(k)k+1,k

...

a(k)i,k 6= 0 . . . ← riga i

...

a(k)n,k . . .

↑colonna k

.

1.1.8 Stabilita della Fattorizzazione L U

Poiche le operazioni aritmetiche che intervengono nell’algoritmo di fatto-rizzazione sono effettuate in aritmetica finita, segue che questi algoritmi,anziche generare i fattori L ed U per A, generano dei fattori non esatti

L = L + δL e U = U + δU.

Posto A + δA = L U , cioe la matrice di cui L U e effettivamente unafattorizzazione, si ha

A + δA = L U = (L + δL)(U + δU) =

= L U + U δL + L δU + δL δU

da cuiδA = U δL + L δU + δL δU

da cui segue che se gli elementi di L ed U sono grandi, gli elementi di δAsono grandi e quindi, gli errori di arrotondamento si amplificano (analisidell’errore all’indietro). Percio diremo che la fattorizzazione A = L U estabile numericamente se gli elementi di L ed U non sono troppo grandirispetto agli elementi di A. In particolare se esistono delle costanti ae b indipendenti dagli elementi e dall’ordine di A tali che |ℓi,j| ≤ a e|ui,j| < b, allora si dice che la fattorizzazione L U e stabile in sensoforte; se le costanti a e b dipendono dall’ordine di A, allora si dice che lafattorizzazione L U e stabile in senso debole.Sebbene si sia visto che per ogni matrice, a meno di una permutazione


di righe, esiste sempre la fattorizzazione L U , in generale questa non estabile. Infatti gli elementi

ℓi,k =a

(k)i,k

a(k)k,k

possono assumere valori grandi e quindi gli elementi di L possono crescereoltre ogni limite.Questo inconveniente puo essere eliminato facendo uno scambio di righead ogni passo dell’algoritmo, cosı da scegliere i pivot a

(k)k,k in modo che

|a(k)k,k| ≥ {|a

(k)i,k |}i=k,...,n.

Questa strategia e indicata con pivoting parziale e limita ad 1 il valoredegli elementi della matrice L, mentre non riesce a limitare gli elementidi U che possono crescere esponenzialmente con l’ordine n di A. Infatti siha che:

max |ui,j| ≤ 2n−1 max |ai,j|.Percio l’algoritmo di Gauss con pivoting parziale genera una fattorizza-zione stabile in senso debole infatti si ha a = 1 e b = 2n−1.

1.1.9 Metodo di Gauss con scambio delle righe eperno massimo

Le considerazioni sopra esposte portano a modificare l’algoritmo di Gaussintroducendo lo scambio delle righe non solo per individuare un pivot nonnullo, ma quello massimo. Questo si realizza al passo k−esimo scegliendocome pivot fra gli a

(k)i,k con i = k, . . . , n il piu grande in valore assoluto,

cioe|a(k)

k,k| ≥ maxi=k+1,...,n

{|a(k)i,k |}.

Quando si usa questo criterio per la scelta dei pivot, il metodo di Gauss sichiama metodo di eliminazione di Gauss con scambio delle righe e perno(o pivot) massimo.

1.2 Condizionamento del Problema Ax = b

In questa sezione si vuole esaminare come perturbazioni sugli elementidella matrice A e sugli elementi del termine noto b influenzano la soluzio-ne x del sistema lineare. Queste perturbazioni sono tipicamente dovuteagli errori di approssimazione quando la matrice A ed il termine noto b

1.2. CONDIZIONAMENTO DEL PROBLEMA AX = B 15

vengono rappresentati con numeri finiti. Per essere in grado di stimare glierrori, dobbiamo introdurre una misura delle dimensioni di un vettore odistanza fra vettori

1.2.1 Richiami sul concetto di norma

Siano x,y ∈ IRn. Una norma vettoriale ‖ · ‖ e una funzione IRn → IRche associa ad un vettore di IRn un valore reale (la lunghezza di quelvettore). Tale funzione, per essere una norma vettoriale deve soddisfarele tre seguenti proprieta:

1. ‖x‖ ≥ 0 ∀x e ‖x‖ = 0 se e solo se x = 0

2. ‖x + y‖ ≤ ‖x‖+ ‖y‖ (disuguaglianza triangolare)

3. ‖αx‖ = |α| ‖x‖ con α ∈ IR.

Le norme vettoriali piu importanti sono:

• ‖x‖1 =∑n

i=1 |xi|

• ‖x‖2 = (∑n

i=1 x2i )

1

2 =√

xTx = (xTx)1

2

• ‖x‖∞ = max1≤i≤n |xi|

Questi sono tutti casi particolari della norma p

‖x‖p =

(

n∑

i=1

|xi|p)

1

p

.

Osservazione 1.1 La ‖ · ‖2 e una generalizzazione ad IRn dell’usuale di-stanza in IR2 o IR3 ed e detta norma Euclidea. La ‖ · ‖∞ e detta normainfinito o norma del massimo.

Con le norme, possiamo introdurre i concetti di distanza e continuitain IRn. Sia x un vettore approssimazione di un vettore x non nullo. Peruna data norma vettoriale ‖ · ‖, si definisceerrore assoluto

‖x− x‖ed errore relativo

‖x− x‖‖x‖ .


Teorema 1.3 (di equivalenza delle norme) Siano ‖ · ‖′ e ‖ · ‖′′ duenorme vettoriali. Allora le due norme sono equivalenti, cioe esistono α eβ ∈ IR con 0 < α ≤ β, tali che per ogni x ∈ IRn e

α ‖ · ‖′′ ≤ ‖ · ‖′ ≤ β ‖ · ‖′′.

Teorema 1.4 Per ogni x ∈ IRn si ha:

• ‖x‖∞ ≤ ‖x‖2 ≤√

n ‖x‖∞• ‖x‖2 ≤ ‖x‖1 ≤

√n ‖x‖2

• ‖x‖∞ ≤ ‖x‖1 ≤ n ‖x‖∞

Sia A una matrice m×n cioe A ∈ IRm×n. Possiamo pensare di definire unanorma di matrici per misurare la dimensione di una matrice e la distanzafra matrici.Una norma matriciale ‖ · ‖ e una funzione IRm×n → IR che associa aduna matrice di IRm×n un valore reale. Tale funzione, per essere una normamatriciale deve soddisfare le tre seguenti proprieta:

1. ‖A‖ ≥ 0 ∀A e ‖A‖ = 0 se e solo se A = 0

2. ‖A + B‖ ≤ ‖A‖+ ‖B‖ (disuguaglianza triangolare)

3. ‖αA‖ = |α| ‖A‖ con α ∈ IR.

Definizione 1.1 (Norma indotta di matrice) Per ogni norma vetto-riale, possiamo definire una corrispondente norma matriciale come

‖A‖ = maxx6=0

‖A x‖‖x‖ .

Osservazione 1.2 Si noti che se A ∈ IRm×n, deve essere x ∈ IRn da cuisi tratta del max di una norma in IRm su una norma in IRn.

Cosı definita, ‖A‖ e una norma di matrice, cioe soddisfa le tre proprietadate.

Risultato 1.1 Se ‖ · ‖ denota una norma vettoriale e la corrispondentenorma matriciale indotta, allora

‖A x‖ ≤ ‖A‖ · ‖x‖

‖A B‖ ≤ ‖A‖ · ‖B‖

1.2. CONDIZIONAMENTO DEL PROBLEMA AX = B 17

Le norme matriciali piu importanti sono:

• ‖A‖1 = max1≤j≤n∑n

i=1 |ai,j|

• ‖A‖∞ = max1≤i≤m∑n

j=1 |ai,j|

• ‖A‖2 =(

λmax(AT A)

)1/2(norma spettrale)

Risultato 1.2 (Relazioni fra le norme matriciali) Per le norme chesono state introdotte valgono le seguenti relazioni, che possono essere di-mostrate dalle relazioni fre le norme vettoriali e la definizione di normaindotta:

• 1√n‖A‖∞ ≤ ‖A‖2 ≤

√n ‖A‖∞

• 1√n‖A‖1 ≤ ‖A‖2 ≤

√n ‖A‖1

• maxi,j |ai,j| ≤ ‖A‖2 ≤ n maxi,j |ai,j|

• ‖A‖2 ≤√

‖A‖1 ‖A‖∞

1.2.2 Errore Inerente

Affrontiamo lo studio dell’errore inerente del problema A x = b conside-rando separatamente eventuali perturbazioni sulla matrice A e sul vettoredei termini noti b.

Introduciamo un vettore di perturbazione δb ∈ IRn sul termine noto;cerchiamo x + δx ∈ IRn soluzione del sistema perturbato

A(x + δx) = b + δb.

poiche e A x = b risulteraAδx = δb

da cuiδx = A−1 δb.

Passando alle norme

‖δx‖ = ‖A−1 δb‖ ≤ ‖A−1‖ · ‖δb‖

inoltre vale‖b‖ = ‖A x‖ ≤ ‖A‖ · ‖x‖.


Allora‖δx‖‖x‖ ≤ ‖A

−1‖ · ‖δb‖‖x‖

≤ ‖A−1‖ · ‖A‖ · ‖δb‖‖b‖ .

Dunque il condizionamento di A x = b dipende dalla costante

K = ‖A‖ · ‖A−1‖

detto numero di condizione di A (rispetto a perturbazioni sul terminenoto).

Introduciamo ora una matrice di perturbazione δA ∈ IRn×n sulla ma-trice dei coefficienti, in modo che A + δA sia ancora invertibile (si puodimostrare che per essere sicuri che A + δA sia invertibile e sufficiente cher = ‖A−1‖ · ‖δA‖ < 1). Cerchiamo x + δx ∈ IRn soluzione del sistemalineare perturbato

(A + δA)(x + δx) = b.

Sviluppando si ha

Ax + δA x + A δx + δA δx = b

δA x + A δx + δA δx = 0

A δx = −δA x− δA δx

δx = −A−1 · δA(x + δx).

Passando alle norme

‖δx‖ = ‖A−1 · δA(x + δx)‖≤ ‖A−1‖ · ‖δA‖ · (‖x‖+ ‖δx‖)

da cui

‖δx‖ ≤ ‖A−1‖ · ‖δA‖ · ‖x‖+ ‖A−1‖ · ‖δA‖ · ‖δx‖

(1− ‖A−1‖ · ‖δA‖)‖δx‖ ≤ ‖A−1‖ · ‖δA‖ · ‖x‖e ricordando che ‖A−1‖ · ‖δA‖ < 1 e ‖A−1‖ · ‖δA‖ = K‖δA‖/‖A‖, siottiene

‖δx‖‖x‖ ≤

K

1− r

‖δA‖‖A‖ ;

dunque ancora il condizionamento di A x = b dipende dalla costante K =‖A‖ · ‖A−1‖ (ora rispetto alla perturbazione sulla matrice dei coefficienti).

1.3. FATTORIZZAZIONE QR 19

Teorema 1.5 Sia A non sigolare e sia r = ‖A−1‖ · ‖δA‖ < 1; allora lamatrice A + δA e non singolare, e

‖(A + δA)−1‖ ≤ ‖A−1‖

1− r.

La soluzione del sistema perturbato

(A + δA)y = b + δb

soddisfa‖y − x‖‖x‖ ≤ K(A)

1− r

(

‖δA‖‖A‖ +

‖δb‖‖b‖

)

.

1.3 Fattorizzazione QR

Data una matrice A di dimensione n×n, esiste sempre una fattorizzazioneQ R con Q matrice ortogonale ed R matrice triangolare superiore:

An×n = Qn×n ·Rn×n

con

QT Q = Iod ancheQT = Q−1

R =

r1,1 . . . r1,n

r2,2

r3,3

. . ....

0 rn,n

.

Osservazione 1.3 Sia Qn×n ortogonale ed a un vettore n × 1, allora‖Qa‖2 = ‖a‖2, infatti

‖Qa‖22 = (Qa)T (Qa) = (aT QT )(Qa) =

= aT (QT Q)a = aT Ia = aTa = ‖a‖22.L’esistenza di una tale fattorizzazione rende facile la determinazio-

ne della soluzione del sistema normale non singolare Ax = b, infatti siavrebbe:

QRx = b

e ponendo Rx = y si potrebbe risolvere il sistema

Qy = b

sfruttando la definizione di Q ortogonale cioe y = QTb ; determinato y sirisolve il sistema

Rx = y

per sostituzione all’indietro (backward-substitution).


1.3.1 Matrici elementari di Householder

Per matrice elementare di Householder si intende una matrice H tale che

H a = ±‖a‖2 e1 =

±‖a‖20...0

.

Limitandoci al caso reale, cioe a vettori ad elementi reali si puo direche le matrici elementari di Householder sono simmetriche (H = HT )ed ortogonali (HT H = I cioe H−1 = HT ).

Assegnato a, l’opportuna trasformazione di Householder si puo costrui-re cosı; si definisce vettore di Householder v come

v =

a1 ± ‖a‖2a2...

an

quindi

H = I − βvvT con β =2

‖v‖22.

Si puo dimostrare che H e ortogonale 1, mentre la simmetria e ovvia dallaforma della H. Inoltre se applicata al vettore a si ha:

H a = (I − βvvT )a= a− βvvTa= a− v infatti βvTa = 1= ∓‖a‖2e1.

Verifichiamo βvTa = 1:

vTa = (a1 ± ‖a‖2, a2, . . . , an)

a1

a2...

an

= a21 ± a1‖a‖2 + a2

2 + . . . + a2n

= ‖a‖22 ± a1‖a‖2= ‖a‖2(‖a‖2 ± a1).

1HT H = (I − βvvT )(I − βvv

T ) = I − 2βvvT + β2

vvTvv

T = I − 2βvvT +

β 2

‖v‖2

2

v‖v‖22v

T = I c.v.d.


inoltre β = 2‖v‖22

e

‖v‖22 = (a1 ± ‖a‖2)2 + a22 + . . . + a2

n)= a2

1 ± 2a1‖a‖2 + ‖a‖22 + a22 + . . . + a2

n

= 2‖a‖22 ± 2a1‖a‖2= 2‖a‖2(‖a‖2 ± a1).

e quindi

βvTa = 1.

Osservazione 1.4 Dato il vettore a = (1, 1, 1)T , si determini la matricedi Householder che azzera la seconda e terza componente.

v =

a1 − ‖a‖2a2

a3

=

1−√

311

β =1

‖a‖2(‖a‖2 − a1)=

1√3(√

3− 1).

H = I − βvvT =

1 0 00 1 00 0 1

− 1√3(√

3− 1)

1−√

311

(1−√

3, 1, 1)

=

1 0 00 1 00 0 1

− 1√3(√

3− 1)

(1−√

3)2 (1−√

3) (1−√

3)

(1−√

3) 1 1

(1−√

3) 1 1

=

1 0 00 1 00 0 1

−

√3− 1√

3− 1√

3− 1√

3− 1√

31√

3(√

3− 1)1√

3(√

3− 1)

− 1√3

1√3(√

3− 1)1√

3(√

3− 1)

=

1 0 00 1 00 0 1

−

1−√

33 −

√3

3 −√

33

−√

33

12 +

√3

612 +

√3

6

−√

33

12 +

√3

612 +

√3

6

=

√3

3

√3

3

√3

3√3

312 −√

36 −1

2 −√

36√

33 −1

2 −√

36

12 −√

36

Ovviamente Ha = (√

3, 0, 0)T .


1.3.2 Metodo di Householder

Procediamo ad illustrare come funziona l’algoritmo per fattorizzare Q Runa matrice An×n = [a1, a2, . . . , an] (il procedimento in oggetto puo essereapplicato anche a matrici rettangolari Am×n con m ≥ n e rango(A) = n).

Si pone A1 := A e si costruisce una matrice elementare di HouseholderH1 tale che

H1 A1 =

‖a(1)1 ‖ × . . . ×0

......

...0 × . . . ×

= A2

quindi una H2 tale che

H2 A2 =

‖a(1)1 ‖ × × . . . ×0 ‖a(2)

2 ‖ × . . . ×... 0

......

0 0 × . . . ×

= A3

e cosı via fino ad una Hn−1 tale che

Hn−1 An−1 =

‖a(1)1 ‖ × × . . . ×0 ‖a(2)

2 ‖ × . . . ×... 0

. . ....

.... . . ×

0 0 . . . ‖a(n−1)n ‖

= R

cosı cheHn−1 Hn−2 . . . H1 A = R.

Ma le Hk sono tutte ortogonali ed il prodotto di matrici ortogonali e ancorauna matrice ortogonale, da cui

A = HT1 HT

2 . . . HTn−1 R = Q R.

Vediamo alcuni dettagli sulla costruzione delle matrici Hk:al primo passo, sia a

(1)1 il vettore formato dagli elementi della prima

colonna di A1 = A e sia

θ1 =

{

+1 se a(1)1,1 ≥ 0

−1 se a(1)1,1 < 0


posto

β1 =1

‖a(1)1 ‖2(‖a(1)

1 ‖2 + |a(1)1,1|)

e

v1 =

θ1(‖a(1)1 ‖2 + |a(1)

1,1|)a

(1)2,1

a(1)3,1...

a(1)n,1

allora la prima matrice elementare di Householder e data da

H1 = I − β1v1vT1 ;

al k−esimo passo, sia a(k)k il vettore di ordine n − k + 1 formato dagli

elementi della k−esima colonna di Ak con indice di riga maggiore e ugualea k, e sia

θk =

+1 se a(k)k,k ≥ 0

−1 se a(k)k,k < 0

posto

βk =1

‖a(k)k ‖2(‖a

(k)k ‖2 + |a(k)

k,k|)e

vk =

0...0

θk(‖a(k)k ‖2 + |a(k)

k,k|)a

(k)k+1,k...

a(k)n,k

la k−esima matrice elementare di Householder e data da

Hk = I − βkvkvTk .

Osservazione 1.5 Se al k−esimo passo si ha a(k)k = 0, si pone Hk = I,

cioe il k−esimo passo non comporta alcuna operazione.


Esempio 1.4 Si calcoli la fattorizzazione Q R della matrice

A1 = A =

72 −144 −144−144 −36 −360−144 −360 450

.

Al primo passo si ha

β1 =1

62208, v1 =

288−144−144

per cui

H1 = I − β1v1vT1 =

1

6=

−2 4 44 4 −24 −2 4

e

A2 =

−216 −216 1080 0 4860 −324 324

.

Al secondo passo si ha

β2 =1

104976, v2 =

0324−324

per cui

H2 = I − β2v2vT2 ==

1 0 00 0 10 1 0

e

R = A3 =

−216 −216 1080 −324 3240 0 −486

.

Inoltre

Q = H1 H2 =1

6

−2 4 44 −2 44 4 −2

.


1.3.3 Implementazione del metodo di Householder

Il metodo di Householder, per risolvere il sistema lineare Ax = b, puo es-sere implementato senza calcolare effettivamente le matrici Hk. Si procedenel seguente modo: si considera la matrice

T1 = [A1|b1] = [A|b]

e si costruiscono β1, v1 ed il vettore di n + 1 componenti

yT1 = vT

1 T1 = (θ1(‖a(1)1 ‖2+|a(1)

1,1|), a(1)2,1, . . . , a

(1)n,1)

a(1)1,1 a

(1)1,2 . . . a

(1)1,n b

(1)1

a(1)2,1 . . .

...... . . .

...

a(1)n,1 . . . a(1)

n,n b(1)n

.

Allora eT2 = H1 T1

= (I − β1v1vT1 ) T1

= T1 − β1v1yT1

Al k−esimo passo si costruisce

yTk = vT

k Tk

e

Tk+1 = Tk − βkvkyTk .

Dopo n− 1 passi si ottiene la matrice

Tn = [An|bn] = [R|bn]

e quindi il sistema Rx = bn con matrice dei coefficienti triangolare supe-riore equivalente al sistema Ax = b.

1.3.4 Costo Computazionale per risolvere Ax = btramite Q R

Si esamina il costo del metodo di soluzione proposto che non comporta ladeterminazione esplicita delle metrici Hk e della matrice ortogonale Q. Alk−esimo passo, per determinare vk e βk, assunte le prime k−1 componenti


di vk nulle e quindi con al piu n− k + 1 componenti non nulle, sara:

vk =

0...0

θk(‖a(k)k ‖2 + |a(k)

k,k|)a

(k)k+1,k...

a(k)n,k

βk =2

‖vk‖22

e servono: n− k + 1 moltiplicazioni ed una radice quadrata per calcolare‖a(k)

k ‖2, una moltiplicazione per ‖vk‖22 ed una divisione per βk per untotale di n− k + 3 molt./div.Per yk servono (n− k + 1)2 + n− k + 1 moltiplicazioni, mentre per vky

Tk

ne servono (n− k + 1)2. In totale

n− k + 3− 2(n− k + 1)2 + (n− k + 1) =

2(n− k)2 + 6(n− k) + 6.

Poiche i passi sono k = 1, . . . , (n− 1) si ha

2(n− 1)n(2n− 1)

6+ 6

(n− 1)n

2+ 4(n− 1) =

2

3(n− 1)n(n + 4) + 4(n− 1) ≈ 2

3n3.

1.3.5 Stabilita Numerica della Fattorizzazione Q R

Procedendo con l’analisi all’indietro, esattamente come gia fatto nel casodella fattorizzazione L U si ha che

δA ∼ R · δQ + Q · δRda cui si cerca una limitazione superiore per gli elementi delle matrici Qed R. Si trova che, essendo Q unitaria, vale

maxi,j|qi,j| ≤ ‖Q‖2 = 1;

inoltre si ha chemax

i,j|ri,j| ≤

√n max

i,j|ai,j|.

Da cui risulta che l’algoritmo di fattorizzazione con matrici elementari diHouseholder (Q R) e stabile in senso debole con estremo

√n, che risulta

comunque piu stabile dell’algoritmo L U essendo√

n≪ 2n−1.

Capitolo 2

Autovalori e Autovettori

Definizione 2.1 data una matrice A ∈ IRm×n, viene definito autovaloredi A un numero λ ∈ IC per cui valga la relazione

A x = λx con x 6= 0. (2.1)

x e detto autovettore corrispondente a λ.

L’insieme degli autovalori di una matrice A costituisce lo spettro diA e l’autovalore massimo in modulo e detto raggio spettrale di A ed eindicato con ρ(A).

La relazione (2.1) data puo essere vista come un sistema lineare omo-geneo

(A− λI) x = 0

che ammette soluzioni non nulle se e solo se

det(A− λI) = 0. (2.2)

Sviluppando la (2.2) risulta

det(A− λI) = p(λ) = a0 + a1λ + . . . + anλn

in cui

an = (−1)n a0 = det(A) an−1 = (−1)n−1tr(A)

dove con tr(A) si indica la somma degli elementi diagonali di A dettatraccia di A.Dalle relazioni che legano i coefficienti e le radici di una equazione algebricarisulta che

n∑

i=1

λi = tr(A) en∏

i=1

λi = det(A).

27

28 CAPITOLO 2. AUTOVALORI E AUTOVETTORI

Il polinomio p(λ) e detto polinomio caratteristico di A e l’equazionep(λ) = 0 e detta equazione caratteristica di A.Per il teorema fondamentale dell’algebra l’equazione caratteristica ha, incampo complesso, n radici. Quindi una matrice di ordine n ha n autova-lori nel campo complesso.Poiche gli autovalori sono soluzioni non nulle del sistema lineare omogeneovisto, un autovettore corrispondente ad un autovalore λ risulta determina-to a meno di una costante moltiplicativa α 6= 0, cioe se x e un autovettoredi A, anche αx e un autovettore di A, corrispondente allo stesso autovalore.

Esempio 2.1 Il polinomio caratteristico della matrice

A =

(

1 33 1

)

si ricava dal determinante

det(A− λI) = det

(

1− λ 33 1− λ

)

= (1− λ)2 − 9 = λ2 − 2λ− 8.

L’equazione caratteristica corrispondente e:

λ2 − 2λ− 8 = 0

ed ha come radici λ1 = −2 e λ2 = 4 che sono gli autovalori della matriceA. L’autovettore corrispondente a λ1 = −2 si calcola risolvendo il sistema

(

3 33 3

)(

x1

x2

)

= 0;

dalla prima equazione si ottiene

x1 + x2 = 0

da cui x1 = −x2 e qualunque vettore

x = α

(

1−1

)

con α 6= 0, e autovettore corrispondente all’autovalore λ1 = −2.L’autovettore corrispondente a λ2 = 4 si determina risolvendo il sistema

(

−3 33 −3

)(

x1

x2

)

= 0;

2.1. PROPRIETA DEGLI AUTOVALORI 29

dalla prima equazione si ottiene

−x1 + x2 = 0

da cui x1 = x2 e qualunque vettore

x = α

(

11

)

con α 6= 0, e autovettore corrispondente all’autovalore λ2 = 4.

Osservazione 2.1 Dato un polinomio si puo costruire una matrice i cuiautovalori sono gli zeri del polinomio; tale matrice si dice companion omatrice di Frobenius. Dato

p(λ) =n∑

i=0

aiλi

la matrice in questione e cosı definita:

F =

0 . . . 0 − a0

an

1. . .

... − a1

an

0 1. . .

......

. . . . . . . . .. . . 0

0 . . . 0 1 −an−1

an

Esempio 2.2 Sia p(λ) = λ2 − 2λ− 8, allora

F =

(

0 81 2

)

verifichiamolo:

det(F − λI) = det

(

−λ −8−1 2− λ

)

= λ(λ− 2)− 8 = λ2 − 2λ− 8.

2.1 Proprieta degli Autovalori

• Gli autovalori di una matrice diagonale o triangolare sono uguali aglielementi diagonali.


• Se λ e un autovalore di una matrice A non singolare e x un auto-vettore corrispondente, allora risulta λ 6= 0 e 1/λ e un autovaloredi A−1 con x autovettore corrispondente. Infatti da Ax = λx si hax = λA−1x e quindi λ 6= 0 e A−1x = 1

λx.

• Se λ e autovalore di una matrice ortogonale A, cioe AT A = AAT = I,allora risulta |λ| = 1. Infatti dalla relazione Ax = λx si ha (Ax)T =(λx)T e quindi xT AT = λxT da cui si ha

xT AT Ax = λxT λx

e poiche A e ortogonale risulta

xTx = λ2xTx e |λ| = 1.

2.2 Proprieta degli Autovettori

Teorema 2.1 Autovettori corrispondenti ad autovalori distinti sono line-armente indipendenti.

Osservazione 2.2 Dal precedente teorema risulta che se una matrice Adi ordine n ha n autovalori tutti distinti, allora A ha n autovettori line-armente indipendenti. Segue che se A non ha autovalori distinti, A puoavere, ma anche non avere, n autovettori linearmente indipendenti.

Definizione 2.2 La molteplicita di un autovalore λ come radice dell’equa-zione caratteristica, e indicata con σ(λ), ed e detta molteplicita algebri-ca di λ. Il massimo numero di autovettori linearmente indipendenti cor-rispondenti a λ e indicato con τ(λ) ed e detta molteplicita geometricadi λ.

Osservazione 2.3 E evidente che

1 ≤ σ(λ) ≤ n e 1 ≤ τ(λ) ≤ n.

Teorema 2.2 Vale la seguente disuguaglianza

τ(λ) ≤ σ(λ).

2.3. SIMILITUDINE FRA MATRICI 31

2.3 Similitudine fra Matrici

Definizione 2.3 Due matrici A e B ∈ IRn×n si dicono simili se esisteuna matrice non singolare S per cui

A = S B S−1.

Osservazione 2.4 La trasformazione che associa la matrice A alla ma-trice B viene detta trasformazione per similitudine; se la matrice S eortogonale, la trasformazione viene detta trasformazione per similitudineortogonale.

Teorema 2.3 Due matrici simili hanno gli stessi autovalori con le stessemolteplicita algebriche e geometriche.

dim. Siano A e B simili, cioe tali che

A = S B S−1.

Si ha che

det(A− λI) = det(S B S−1 − λS S−1)= det(S(B − λI)S−1)= det(S) det(B − λI) det(S−1)= det(B − λI)

per cui le due matrici hanno lo stesso polinomio caratteristico e quindihanno gli stessi autovalori con le stesse molteplicita algebriche.Se x e autovettore di A corrispondente all’autovalore λ, risulta

S B S−1x = λx

e quindiB S−1x = λS−1x

percio il vettore y = S−1x e autovettore di B corrispondente a λ.Inoltre, essendo S−1 non singolare, se xi, i = 1, . . . , τ(λ) sono autovettorilinearmente indipendenti di A, anche yi = S−1xi, i = 1, . . . , τ(λ) sonolinearmente indipendenti.

Definizione 2.4 Una matrice A simile ad una matrice diagonale D sidice diagonalizzabile.

Teorema 2.4 Una matrice A di ordine n e diagonalizzabile se e solo seha n autovettori linearmente indipendenti. Inoltre le colonne della matriceS, per cui S−1 A S e diagonale, sono gli autovettori di A.


dim. Siano x1,x2, . . . ,xn n autovettori linearmente indipendenti di Acorrispondenti agli autovalori λ1, λ2, . . . , λn. Siano D la matrice diagonaleavente λi come i−esimo elemento diagonale, ed S la matrice la cui i−esimacolonna e uguale a xi. Dalla relazione

Axi = λixi i = 1, . . . n

si ha cheA S = S D.

Essendo S non singolare, perche formata da colonne linearmente indipen-denti, esiste S−1; quindi si ha

A = S D S−1.

Viceversa, sia A = S D S−1, con D matrice diagonale con gli autovaloridi A come elementi diagonali. Allora risulta

A S = S D.

Indicando con s1, s2, . . . , sn le colonne di S, si ha

A[s1, s2, . . . , sn] = [λ1s1, λ2s2, . . . , λnsn].

Percio le colonne di S sono n autovettori di A che risultano linearmenteindipendenti.

Esempio 2.3 Sia

A =

1 1 01 1 10 1 1

.

Questa ha come autovalori

λ1 = 1−√

2 λ2 = 1 λ3 = 1 +√

2

e i corrispondenti autovettori sono

x1 =

1

−√

21

x2 =

10−1

x3 =

1√21

.

Allora A e diagonalizzabile, cioe esiste una matrice S, non singolare, taleche

A = S D S−1.

2.3. SIMILITUDINE FRA MATRICI 33

Infatti per quanto visto sara:

S =

1 1 0

−√

2 0√

20 1 1

da cui risulta

S−1 =

1/4 −√

2/4 1/41/2 0 −1/2

1/4√

2/4 1/4

e

A = S

1−√

2 0 00 1 0

0 0 1 +√

2

S−1.

Osservazione 2.5 Fra le trasformazioni per similitudine che associanoalla matrice B la matrice A = S B S−1, hanno particolare importanzaquelle per cui S e ortogonale, cioe ST S = S ST = I. Il teorema che se-gue mostra come sia possibile, mediante una trasformazione per similitudi-ne ortogonale, ricondurre una qualsiasi matrice ad una forma triangolaresuperiore.

Teorema 2.5 (Forma Normale di Schur) Sia A ∈ ICn e siano λ1, . . . , λn

i suoi autovalori. Allora esiste una matrice ortogonale/unitaria Q e unamatrice triangolare superiore T i cui elementi diagonali sono i λi, tali che

A = Q T QT .

se la matrice A ha elementi reali esiste la forma normale reale di Schur.

Teorema 2.6 (Forma Normale Reale di Schur) se A ∈ IRn×n, esisteuna matrice ortogonale Q ∈ IRn×n e una matrice T ∈ IRn×n triangolaresuperiore a blocchi della forma

T =

R1,1 R1,2 . . . R1,m

R2,2...

. . .

0 Rm,m

.

dove i blocchi Rj,j per j = 1, . . . ,m hanno ordine 1 o 2. Se λj e autovalorereale di A, allora Rj,j ha ordine 1 e coincide con λj, se invece e complesso,allora il blocco Rj,j ha ordine 2 ed ha come autovalori λj e λj (il complessoconiugato). La somma delle dimensioni dei blocchi Rj,j j = 1, . . . ,m e pariad n.


Osservazione 2.6 Una classe particolarmente importante di matrici equella delle matrici normali, cioe tali che

A AH = AH A

dove H sta per Hermitiana ed equivale ad una matrice simmetrica in cam-po complesso. Questa classe e importante perche comprende tutte e sole lematrici diagonalizzabili con trasformazioni per similitudine unitarie. Valeinfatti il seguente teorema.

Teorema 2.7 Una matrice A ∈ ICn×n e normale (cioe A AH = AH A) see solo se esiste una matrice unitaria Q tale che

A = Q

λ1

λ2

. . .

λn

QT

in cui λ1, λ2, . . . , λn sono gli autovalori di A. La matrice Q ha per colonnegli autovettori di A, che quindi sono a due a due ortonormali.

Capitolo 3

Sistemi lineari: metodiiterativi

Per risolvere un sistema lineare Ax = b, oltre ai metodi diretti, si possonoutilizzare anche i metodi iterativi, che risultano particolarmente conve-nienti se la matrice e sparsa, cioe se il numero degli elementi non nullidi A e dell’ordine della dimensione della matrice. Infatti quando si usaun metodo diretto puo accadere che nelle matrici intermedie vengano ge-nerati molti elementi diversi da zero in corrispondenza di elementi nullidella matrice iniziale. Poiche i metodi diretti non sfruttano la sparsita del-la matrice, soprattutto poi se A e anche di grandi dimensioni, puo essereconveniente utilizzare un metodo iterativo. Esistono comunque dei casi neiquali la matrice A e sparsa, ma e piu conveniente usare dei metodi direttiche sfruttano specifiche proprieta di struttura della matrice. Cominciamocon alcuni richiami su successioni convergenti di vettori e matrici.

Definizione 3.1 Una successione {x(k)} di vettori di IRn si dice conver-gente al vettore x∗ ∈ IRn se esiste una norma vettoriale per cui

limk→∞‖x(k) − x∗‖ = 0;

in tal caso si ponelim

k→∞x(k) = x∗.

Per il teorema di equivalenza delle norme, la definizione appena data nondipende da una particolare norma. La condizione di convergenza data sitraduce in una condizione di convergenza delle successioni formate dallesingole componenti. Infatti, considerando la norma ∞ o del max

|x(k)i − x∗

i | ≤ ‖x(k) − x∗‖∞ i = 1, . . . , n

35

36 CAPITOLO 3. SISTEMI LINEARI: METODI ITERATIVI

e quindilimk→∞|x(k)

i − x∗i | = 0

da cuilim

k→∞x

(k)i = x∗

i .

Viceversa se vale quest’ultima e ovviamente verificata la condizione diconvergenza data nella definizione, per la norma ∞.

Per le successioni di matrici {A(k)} si puo dare una definizione di con-vergenza analoga a quella data per vettori. Il seguente teorema e di fon-damentale importanza nello studio dei metodi iterativi per la soluzione disistemi lineari.

Teorema 3.1 Sia A ∈ IRn×n, allora

limk→∞

Ak = 0 se e solo se ρ(A) < 1

con ρ(A) il raggio spettrale della matrice A.

Teorema 3.2 Sia A ∈ IRn×n, allora

det(I − A) 6= 0

e

limk→∞

k∑

i=0

Ai = (I − A)−1 se e solo se ρ(P ) < 1.

Osservazione 3.1 Come per le serie numeriche, si usa scrivere

∞∑

i=0

Ai = (I − A)−1.

3.1 Decomposizione della matrice

Sia A ∈ IRn×n una matrice non singolare e si consideri la decomposizionedi A nella forma

A = M − N

dove M e una matrice non singolare. Sostituendo tale decomposizione nelsistema Ax = b si ha

(M − N)x = b

Mx − Nx = b

3.1. DECOMPOSIZIONE DELLA MATRICE 37

ed essendo M non singolare

x = M−1 Nx + M−1b.

Posto P := M−1 N e q := M−1b si ottiene il seguente sistema

x = Px + q

equivalente ad Ax = b.Dato un vettore iniziale x(0), si considera la successione x(1),x(2), . . . cosıdefinita

x(k) = Px(k−1) + q k = 1, 2, . . . (3.1)

Se la successione {x(k)} e convergente, cioe

limk→∞

x(k) = x∗

allora passando al limite nella (3.1) risulta

x∗ = Px∗ + q (3.2)

cioe x∗ e la soluzione del sistema x = Px+q e quindi del sistema Ax = b.La relazione (3.1) individua un metodo iterativo in cui, partendo da unvettore iniziale x(0), la soluzione viene approssimata utilizzando una suc-cessione x(k) di vettori. La matrice P si dice matrice di iterazione delmetodo.

Osservazione 3.2 Al variare del vettore iniziale x(0) si ottengono dalla(3.1) diverse successioni {x(k)}, alcune delle quali possono essere conver-genti ed altre no. Un metodo iterativo e detto convergente se, qualunquesia il vettore iniziale x(0), la successione {x(k)} e convergente.

Teorema 3.3 Il metodo iterativo (3.1) risulta convergente se e solo seρ(A) < 1.

dim. Sottraendo dalla (3.2) la (3.1) si ha

x∗ − x(k) = P (x∗ − x(k−1)) k = 1, 2, . . .

Indicato cone(k) = x∗ − x(k)

il vettore errore alla k−esima iterazione, si ha

e(k) = Pe(k−1) k = 1, 2, . . .


e quindie(k) = Pe(k−1) = P 2e(k−2) = . . . = P ke(0).

(Condizione sufficiente) Se ρ(P ) < 1, per il teorema visto, risulta

limk→∞

P k = 0

e segue che per ogni vettore e(0) si ha

limk→∞

e(k) = 0.

(Condizione necessaria) Se il metodo e convergente il limk→∞ e(k) = 0vale per ogni x(0), e in particolare deve valere se x(0) e tale che il vettotee(0) = x∗ − x(0) e un autovettore di P corrispondente ad un autovalore λdi modulo massimo, cioe |λ| = ρ(P ).In questo caso risulta

Pe(0) = λe(0)

e quindie(k) = P ke(0) = λke(0).

Segue chelim

k→∞[ρ(P )]k = 0

e quindi ρ(P ) =< 1.La condizione ρ(P ) < 1, necessaria e sufficiente per la convergenza

della (3.1), non e in generale di agevole verifica. Conviene utilizzare,quando e possibile, delle condizioni sufficienti di convergenza di piu facileverifica. Una tale condizione e data dal seguente teorema.

Teorema 3.4 Se esiste una norma matriciale indotta ‖·‖, per cui ‖P‖ <1, il metodo iterativo (3.1) e convergente.

dim. La tesi segue dal teorema di convergenza visto e dalla proprieta che

ρ(P ) < ‖P‖.

Osservazione 3.3 Il raggio spettrale ρ(P ) di una matrice viene a volteanche definito come l’estremo inferiore di tutte le norme di P .

Osservazione 3.4 Poiche il determinante di una matrice e uguale al pro-dotto degli autovalori, se |det(P )| ≥ 1, almeno uno degli autovalori di Pe in modulo magiore o uguale ad 1 e quindi il metodo (3.1) non e conver-gente.

3.2. CONTROLLO DELLA CONVERGENZA 39

Poiche la traccia di una matrice e uguale alla somma degli autovalori, se|tr(P )| > n, almeno uno degli autovalori di P e in modulo maggiore ouguale ad 1 e il metodo (3.1) non e convergente.Quindi le condizioni:

|det(P )| < 1 e |tr(P )| < n

sono condizioni necessarie affinche il metodo (3.1) sia convergente.

3.2 Controllo della Convergenza

Se e(k−1) 6= 0, la quantita

‖e(k)‖‖e(k−1)‖

esprime la riduzione dell’errore al k−esimo passo.la media geometrica delle riduzioni dell’errore sui primi k passi

σk =

(

‖e(1)‖‖e(0)‖

· ‖e(2)‖

‖e(1)‖· · · ‖e

(k)‖‖e(k−1)‖

)

1

k

=

(

‖e(k)‖‖e(0)‖

)

1

k

esprime la riduzione media per passo dell’errore relativo ai primi kpassi. Dalla

‖e(k)‖ ≤ ‖P k‖ · ‖e(0)‖

risulta

σk ≤ (‖P k‖) 1

k .

La quantita che si ottiene facendo tendere k all’infinito esprime la ridu-zione asintotica media per passo e, come risulta dal seguente teorema,e indipendente dalla particolare norma utilizzata.

Teorema 3.5 Sia A ∈ IRn×n e sia ‖ · ‖ una qualunque norma matricialeindotta. Allora

limk→∞

(‖P k‖) 1

k = ρ(P ).

Osservazione 3.5 La quantita ρ(P ), indipendente dalla norma utilizzatae dall’indice di iterazione k, viene quindi assunta come misura dellavelocita di convergenza del metodo (3.1).


3.3 Test di Arresto

Poiche con un metodo iterativo non e possibile calcolare la soluzione conun numero finito di iterazioni, occorre individuare dei criteri per l’arrestodel procedimento.I criteri piu comunemente usati, fissata una tolleranza Tol, che tiene contoanche della precisione utilizzata nei calcoli, sono i seguenti

‖x(k) − x(k−1)‖ ≤ Tol

oppure, se x(k) 6= 0

‖x(k) − x(k−1)‖‖x(k)‖

≤ Tol.

Si noti che queste due condizioni non garantiscono che la soluzione siastata approssimata con la precisione Tol. Infatti da

e(k) = Pe(k−1)

ex(k) − x(k−1) = (x∗ − x(k−1))− (x∗ − x(k))

= e(k−1) − e(k)

= (I − P )e(k−1)

e passando alle norme, se ‖P‖ < 1, per il teorema seguente si ha

‖e(k−1)‖ ≤ ‖(I − P )−1‖ · ‖x(k) − x(k−1)‖≤ ‖x(k) − x(k−1)‖

1− ‖P‖ .

Per cui puo accadere che ‖e(k−1)‖ sia elevata anche se le condizioni diarresto date sono verificate.

Teorema 3.6 Sia A ∈ IRn×n e sia ‖ · ‖ una qualunque norma matricialeindotta tale che ‖A‖ < 1. Allora la matrice I ± A e non singolare e valela disuguaglianza

‖(I ± A)−1‖ ≤ 1

1− ‖A‖ .

• In un programma che implementa un tale metodo iterativo deveessere comunque previsto un controllo per interrompere l’esecuzionequando il numero delle iterazioni diventa troppo elevato.

3.4. METODI DI JACOBI E GAUSS-SEIDEL 41

• Puo anche accadere che un metodo iterativo la cui matrice di ite-razione P e tale che ρ(P ) < 1, per gli effetti indotti dagli errori diarrotondamento non converga in pratica. Questo accade, in parti-colare, quando la matrice A e fortemente mal condizionata e ρ(P ) emolto vicino ad 1.

• Si deve rilevare che un metodo iterativo, rispetto ad un metodo diret-to e in generale meno sensibile alla propagazione degli errori. Infattiil vettore x(k) puo essere considerato come il vettore generato conuna sola iterazione a partire dal vettore iniziale x(k−1), e quindi risul-ta affetto dagli errori di arrotondamento generati dalla sola ultimaiterazione.

• In un metodo iterativo, ad ogni iterazione, il costo computazionale eprincipalmente determinato dalle operazioni di moltiplicazione dellamatrice P per un vettore, che richiede n2 operazioni moltiplicative.Se A e sparsa, il numero di moltiplicazioni e dell’ordine di n. Inquesto caso i metodi iterativi possono risultare competitivi con quellidiretti.

3.4 Metodi di Jacobi e Gauss-Seidel

Tra i metodi iterativi individuati da una particolare scelta della decompo-sizione A = M −N sono particolarmente importanti il metodo di Jacobie il metodo di Gauss-Seidel, per i quali e possibile dare delle condizionisufficienti di convergenza. Si consideri la decomposizione della matrice A

A = D − L− U

con

D =

a1,1

a2,2

. . .

an,n

,

L =

0 . . . 0

−a2,1. . .

......

. . .

−an,1 . . . −an,n−1 0

, U =

0 −a1,2 . . . −a1,n...

. . . . . .

−an−1,n

0 . . . 0

.

ScegliendoM = D N = L + U,


si ottiene il metodo di Jacobi, mentre scegliendo

M = D − L N = U,

si ottiene il metodo di Gauss-Seidel.Per queste decomposizioni risulta det(M) 6= 0 se e solo se tutti gli elementidiagonali di A sono non nulli.Indicando con J la matrice di iterazione del metodo di Jacobi, dalla P =M−1N si ha

J = D−1(L + U)

per cui la (3.1) diviene

x(k) = Jx(k−1) + D−1b = D−1[b + (L + U)x(k−1)]

e in termini di componenti

x(k)i =

1

ai,i

bi −n∑

j=1 j 6=i

ai,jx(k−1)j

i = 1, . . . , n.

Il metodo di Jacobi e detto anche metodo degli spostamenti simultanei,in quanto le componenti del vettore x(k) sostituiscono simultaneamente,al termine dell’iterazione, le componenti di x(k−1).

Indicando con G la matrice di iterazione di Gauss-Seidel, dalla P =M−1N si ha

G = (D − L)−1U

per cui la (3.1) diviene

x(k) = Gx(k−1) + (D − L)−1b.

Per descrivere questo metodo in termini di componenti, conviene primatrasformarla nel modo seguente:

(D − L)x(k) = Ux(k−1) + b

Dx(k) = Lx(k) + Ux(k−1) + b

x(k) = D−1(

Lx(k) + Ux(k−1) + b)

ottenendo quindi

x(k)i =

1

ai,i

bi −i−1∑

j=1

ai,jx(k)j −

n∑

j=i+1

ai,jx(k−1)j

i = 1, . . . , n.

3.4. METODI DI JACOBI E GAUSS-SEIDEL 43

Confrontando questa con quella di Jacobi, risulta che per calcolare le com-ponenti del vettore x(k) sono utilizzate componenti gia calcolate dello stes-so vettore. Per questo motivo il metodo di Gauss-Seidel e detto anchemetodo degli spostamenti successivi.

Nell’implementare Jacobi e necessario disporre, contemporaneamentedi entrambe i vettori x(k−1) e x(k), mentre per Gauss-Seidel e sufficientedisporre di un solo vettore.

Osservazione 3.6 In molte applicazioni il metodo di Gauss-Seidel, cheutilizza immediatamente i valori calcolati nella iterazione corrente, risultapiu veloce del metodo di Jacobi. Pero esistono casi in cui risulta non soloche il metodo di Jacobi sia piu veloce di Gauss-Seidel, ma anche che Jacobisia convergente e Gauss-seidel no.

Per i metodi di Jacobi e Gauss-Seidel si possono ricavare delle condizionisufficienti di convergenza di facile verifica sul sistema lineare.

Teorema 3.7 Se la matrice A e a predominanza diagonale in senso stret-to, allora i metodi di Jacobi e di Gauss-Seidel sono convergenti.

Osservazione 3.7 Una matrice A ∈ IRn×n si dice a predominanza diago-nale in senso stretto, se vale

|ai,j| >n∑

j=1 j 6=i

|ai,j| i = 1, . . . , n per righe

|ai,j| >n∑

i=1 i6=j

|ai,j| j = 1, . . . , n per colonne

Esempio 3.1 Applichiamo il metodo di Jacobi al sistema lineare

Ax = b con A =

20 2 −12 13 −21 1 1

e b =

25302

.

La decomposizione di Jacobi e

M =

20 0 00 13 00 0 1

, N =

0 −2 1−2 0 2−1 −1 0

La matrice di iterazione J ed il vettore q sono

J =

0 −1/10 1/20−2/13 0 2/13−1 −1 0

, q =

5/430/13

2


Risulta ρ(J) = 0.44896 < 1, da cui il metodo converge e converge allasoluzione x = (1, 2,−1)T . Sia x(0) = (0, 0, 0)T , allora:

x(1) = (1.25, 2.3 . . . , 2)T ‖e(1)‖2 ≈ 3.026x(2) = (1.12 . . . , 2.42 . . . ,−1.55 . . .)T ‖e(2)‖2 ≈ 0.71x(3) = (0.93 . . . , 1.89 . . . ,−1.54 . . .)T ‖e(3)‖2 ≈ 0.55x(4) = (0.98 . . . , 1.92 . . . ,−0.82 . . .)T ‖e(4)‖2 ≈ 0.19x(5) = (1.01 . . . , 2.02 . . . ,−0.91 . . .)T ‖e(5)‖2 ≈ 0.09x(6) = (1.00 . . . , 2.01 . . . ,−1.04 . . .)T ‖e(6)‖2 ≈ 0.046

Esempio 3.2 Applichiamo il metodo di Gauss-Seidel al sistema lineare

Ax = b con A =

20 2 −12 13 −21 1 1

e b =

25302

.

La decomposizione di Gauss-Seidel e

M =

20 0 02 13 01 1 1

, N =

0 −2 10 0 20 0 0

.

La matrice di iterazione G ed il vettore q sono

G =

0 −1/10 1/200 1/65 9/1300 11/130 −51/260

, q =

5/455/26−71/52

Risulta ρ(G) = 0.243858 < 1, da cui il metodo converge e converge allasoluzione x = (1, 2,−1)T . Sia x(0) = (0, 0, 0)T , allora:

x(1) = (1.25, 2.11 . . . ,−1.36)T ‖e(1)‖2 ≈ 0.45x(2) = (0.97 . . . , 2.05 . . . ,−0.91 . . .)T ‖e(2)‖2 ≈ 0.101

x(3) = (0.998 . . . , 2.08 . . . ,−1.011 . . .)T ‖e(3)‖2 ≈ 0.084x(4) = (0.991 . . . , 2.07 . . . ,−0.99 . . .)T ‖e(4)‖2 ≈ 0.078

Capitolo 4

Metodi per Autovalori eAutovettori

Tra i diversi metodi numerici esistenti per calcolare gli autovalori e gli au-tovettori di una matrice, alcuni sono di tipo generale e sono applicabili amatrici dense e senza struttura, altri utilizzano in modo specifico eventua-li proprieta di struttura o sparsita della matrice, permettendo di trattareproblemi anche con dimensioni molto grandi. Alcuni metodi possono esse-re utilizzati per calcolare tutti gli autovalori e autovettori di una matrice,altri invece servono per calcolare solo alcuni autovalori, per esempio quelliche si trovano all’estremita dello spettro, ed i corrispondenti autovettori,come e richiesto in molte applicazioni. Tutti i metodi possono essere divisiin due classi:

1. Metodi in cui il calcolo viene effettuato in due fasi

• riduzione con metodi diretti della matrice A in una matricesimile B, di cui sia piu agevole calcolare gli autovalori;

• calcolo degli autovalori di B con un metodo iterativo.

Questi metodi si applicano in generale a problemi di piccole dimen-sioni, per i quali tutti i dati su cui si opera possono essere contenutinella memoria centrale del calcolatore.

2. Metodi completamente iterativi, che richiedono ad ogni passo la mol-tiplicazione di una matrice per un vettore, o la risoluzione di un si-stema lineare. Questi metodi si applicano in generale a problemi digrandi dimensioni, anche nel caso in cui non sia possibile conteneretutti i dati nella memoria del calcolatore.

45

46 CAPITOLO 4. METODI PER AUTOVALORI E AUTOVETTORI

Nei metodi della prima classe, per la riduzione della matrice A nella ma-trice B, si utilizzano metodi diretti analoghi a quelli descritti per la fatto-rizzazione delle matrici. Nel caso piu generale la matrice B che si ottienee tale che

bi,j = 0 per i > j + 1 i, j = 1, . . . , n.

Una matrice B con questa proprieta e detta essere in forma di Hessen-berg superiore. La trasformazione per similitudine della matrice A nellamatrice B e fatta per passi successivi

Ak+1 = Sk Ak S−1k k = 1, . . . ,m− 1

doveA1 := A e alla fine B := Am

per cui posto S := Sm−1 Sm−2 . . . S1, risulta

B = S A S−1

e se x e autovettore di B, Sx e autovettore di A.Le matrici Sk sono di solito matrici elementari di Gauss o di Householder.Se Sk e una matrice di Householder, risulta

‖Sk‖2 · ‖S−1k ‖2 = 1

I metodi iterativi per il calcolo degli autovalori di B potrebbero essereanche applicati direttamente alla matrice A. Trasformando pero la matriceA nella matrice B, si abbassa notevolmente la complessita computazionale(per esempio per il metodo QR che vedremo si passa da O(n3) a O(n2)).Per il calcolo degli autovalori della matrice B, due sono le tecniche piuusate:

• se sono richiesti solo pochi autovalori rispetto alla dimensione dellamatrice (non piu del 25%), conviene usare un metodo iterativo checalcoli un singolo autovalore per volta, come ad esempio un meto-do di iterazione funzionale applicato all’equazione caratteristica o ilmetodo delle potenze inverse.

• se sono richiesti tutti o molti degli autovalori, il metodo migliore ein generale il metodo QR che vedremo.

Osservazione 4.1 Se la matrice A e hermitiana, e la trasformazione vie-ne eseguita con matrici unitarie, la matrice B risulta hermitiana tridia-gonale.

4.1. RIDUZIONE DI UNA MATRICE IN FORMA DI HESSENBERG47

4.1 Riduzione di una matrice in forma di

Hessenberg

Se si applicano ad una matrice An×n matrici elementari di Householdercon l’obiettivo di azzerare gli elementi della colonna j−esima di indicii = j +2, . . . , n si ottiene una matrice di Hessenberg superiore ad un costocomputazionale di 5

3n3 operazioni moltiplicative.

Esempio 4.1 Si consideri la matrice A ∈ IR4×4

A =

4 3 2 11 4 3 21 1 4 31 1 1 4

Applicando il metodo di Householder, al primo passo si ottiene

v1 = (0.0, 2.732051, 1.0, 1.0)T β1 = 0.2113248

e quindi H1 = I − β1v1vT1 e A2 = H1AHT

1

A2 =

4 −3.464098 −0.366024 −1.366024−1.73205 7.666641 −0.5446615 −1.12209

0 0.08931351 1.877991 1.0326920 1.244013 −0.6993591 2.455341

.

Al secondo passo si ottiene

v2 = (0.0, 0.0, 1.336528, 1.244013)T β2 = 0.5999027

e quindi H2 = I − β2v2vT2 e A3 = H2A2H

T2

A3 =

4 −3.464098 1.388726 0.2672625−1.73205 7.666641 1.158131 0.4629154

0 −1.247213 2.476189 −0.74230770 0 0.9897442 1.857142

che e in forma di Hessenberg superiore.

4.2 Metodo Q R per il calcolo degli autova-

lori

Il metodo QR e il metodo piu usato per calcolare tutti gli autovalori diuna matrice, in quanto il piu efficiente.


Il metodo e molto complicato, sia come descrizione che come implemen-tazione, anche se il principio su cui si basa e semplice. Il metodo richiedetutta una serie di accorgimenti, senza i quali non potrebbe essere efficien-te: riduzione preliminare della matrice in forma di Hessenberg superiore(o triangolare nel caso di matrici hermitiana) per ridurre il costo compu-tazionale ad ogni iterazione, utilizzazione di una tecnica di traslazione peraumentare la velocita di convergenza; riduzione dell’ordine della matricequando un autovalore e stato approssimato con sufficiente precisione, percalcolare un altro autovalore.Il metodo QR, che e stato descritto da Francis nel 1961, utilizza la fatto-rizzazione QR di una matrice; esso deriva da un precedente metodo, dettometodo LR, proposto da Rutishauser nel 1958, che utilizza la fattorizza-zione LU di una matrice.

4.2.1 Algoritmo di base

Nel metodo QR viene generata una successione {Ak} di matrici nel modoseguente: posto

A1 := A

per k = 1, 2, . . . si calcola una fattorizzazione QR di Ak

Ak = QkRk

dove Qk e unitaria e Rk e triangolare superiore, e si definisce la matriceAk+1 per mezzo della relazione

Ak+1 := RkQk.

per le precedenti risulta che

Ak+1 = QTk AkQk

e quindi le matrici della successione {Ak} sono tutte simili fra loro. Sottoopportune ipotesi la successione converge ad una matrice triangolare su-periore (diagonale se A e hermitiana) che ha come elementi diagonali gliautovalori di A.

4.2.2 Risultati di convergenza

Il seguente teorema di convergenza, viene dato in ipotesi piuttosto restrit-tive in cui e piu facile la sua dimostrazione, ma la sua validita puo essereprovata anche in ipotesi piu deboli.

4.2. METODO Q R PER IL CALCOLO DEGLI AUTOVALORI 49

Teorema 4.1 Sia A ∈ ICn×n tale che i suoi autovalori λi, i = 1, . . . , nabbiano moduli tutti distinti, cioe

|λ1| > |λ2| > . . . > |λn| > 0.

Indicata con X la matrice degli autovettori di A, tale che

A = XDX−1,

in cui D e la matrice diagonale degli autovalori, si supponga che la matriceX−1 ammetta la fattorizzazione LU . Allora esistono delle matrici Sk taliche

limk→∞

SHk RkSk−1 = lim

k→∞SH

k−1AkSk−1 = T

e

limk→∞

SHk−1QkSk = I

dove T e triangolare superiore con gli elementi diagonali uguali a λ1, λ2, . . . , λn.Quindi gli elementi diagonali di Ak tendono agli autovalori di A.

Se X−1 non ammette fattorizzazione LU , si puo dimostrare che il metodoQR e ancora convergente. In questo caso gli elementi diagonali di T coin-cidono ancora con i λi, ma non sono piu in ordine di modulo decrescente.Se l’ipotesi che tutti gli autovalori abbiano modulo distinto, non e verifi-cata, la successione formata dagli elementi diagonali di Ak non converge.Questa ipotesi e troppo restrittiva, e non consente di utilizzare il metodoQR in casi particolarmente importanti nelle applicazioni, come quelli incui la matrice A ha elementi reali e autovalori non reali. Pero anche inquesto caso il metodo QR puo essere applicato con opportune varianti.Sia ad esempio

|λ1| > . . . > |λr| = |λr+1| > . . . > |λn| > 0

dove λr e λr+1 sono due numeri complessi coniugati, oppure due numerireali. Allora la successione delle {Ak} o meglio degli elementi diagonalinon converge agli autovalori di stesso modulo, ma gli autovalori dei blocchidiagonali convergono a λr e λr+1. Come detto, gli elementi diagonali diAk di indice diverso da r ed r + 1 invece convergono agli altri autovalori.Situazioni analoghe si presentano quando la matrice A ha piu autovaloridi modulo uguale e in questo caso il metodo QR genera matrici Rk construttura triangolare a blocchi, in cui gli autovalori dei blocchi diagonaliconvergono ad autovalori di A.


4.2.3 Costo computazionale e stabilita

Il metodo QR applicato a una matrice di ordine n ha ad ogni passo uncosto computazionale dell’ordine di n3 operazioni moltiplicative (per calco-lare la fattorizzazione Ak = QkRk e per moltiplicare la matrice triangolareRk per le matrici elementari di Householder).Per abbassare il costo computazionale globale conviene prima trasformarela matrice A in forma di Hessenberg superiore. Questa trasformazioneviene eseguita una sola volta perche il metodo QR, applicato a matrici informa di Hessenberg superiore produce matrici Ak in forma di Hessenbergsuperiore.Infatti se Ak e in forma di Hessenberg superiore, la matrice Qk e data dalprodotto di n − 1 matrici elementari di Hoseholder che sono in forma diHessenberg superiore e quindi la matrice Ak+1, prodotto di una matricetriangolare superiore Rk per una matrice Qk in forma di Hessenberg su-periore, risulta ancora in forma di Hessenberg superiore.Se la matrice A e hermitiana, la matrice in forma di Hessenberg superiore,ottenuta applicando ad A il metodo di Householder, e ancora hermitiana,e quindi risulta tridiagonale. Inoltre anche tutte le matrici Ak generatedal metodo QR sono hermitiane e quindi tridiagonali.Il metodo QR applicato ad una matrice A in forma di Hessenberg supe-riore ha ad ogni passo un costo di 2n2 operazioni moltiplicative (che e ilcosto computazionale per calcolare la fattorizzazone Ak = QkRk).Per quel che riguarda la stablita si puo dimostrare che il metodo QR per ilcalcolo degli autovalori gode delle stesse proprieta di stabilita di cui godela fattorizzazione QR di una matrice.

4.3 Condizionamento del calcolo degli au-

tovalori

Una matriice puo avere alcuni autovalori molto sensibili alle perturbazionidei suoi elementi mentre altri insensibili. E conveniente avere un numeroche definisce il condizionamento di una matrice rispetto al problema degliautovalori; esso sara il numero di condizione spettrale.E evidente che tale numero puo fornire una informazione insoddisfacente.Infatti se qualche autovalore e molto sensibile allora il numero di condi-zione spettrale e necessariamente grande, anche se altri autovalori sonoinsensibili. Si puo ovviare a questo inconveniente introducendo numeridi condizione che evidenziano la sensibilita dei singoli autovalori e che si

4.3. CONDIZIONAMENTO DEL CALCOLO DEGLI AUTOVALORI51

chiamano numeri di condizione della matrice rispetto al problemadegli autovalori.

Definizione 4.1 Si dice norma assoluta ‖ · ‖ una norma matriciale in-dotta che verifichi le seguenti proprieta

‖D‖ = maxi=1,...,n

|di,i|

per ogni matrice diagonale D ∈ ICn×n.

Osservazione 4.2 Le norme ‖ · ‖1, ‖ · ‖2 e ‖ · ‖∞ sono assolute.

Teorema 4.2 (di Bauer-Fike) Sia A ∈ ICn×n una matrice diagonalizza-bile, cioe tale che

A = TDT−1

con D diagonale e T non singolare. Se δA ∈ ICn×n e η e un autovalore diA + δA, allora esiste almeno un autovalore λ di A tale che

|λ− η| ≤ µ(T )‖δA‖

dove µ(T ) = ‖T‖ · ‖T−1‖ per una norma assoluta ‖ · ‖.

Introduciamo un numero di condizione che dia una misura dell’effetto dellaperturbazione su un solo autovalore di A.

Teorema 4.3 Sia A ∈ ICn×n, λ un autovalore di A di molteplicita alge-brica 1, x,y ∈ ICn×n, ‖x‖2 = ‖y‖2 = 1, tali che

Ax = λx

yHA = λyH .

Allora e yHx 6= 0 ed inoltre esiste nel piano complesso un intorno V dellozero e una funzione λ(ǫ) : V → IC, analitica, tale che

1. λ(ǫ) e autovalore con molteplicita 1 di A + ǫF , F ∈ ICn×n

2. λ(0) = λ

3. λ′(0) = yHFx

yHx

4. a meno di termini di ordine superiore in ǫ e

λ(ǫ)− λ = ǫyHFx

yHx.


da questo risultato si ha che la variazione nell’autovalore dovuta alla per-turbazione ǫF di A e proporzionale ad ǫ.Inoltre il condizionamento del problema dipende dalla quantita

∣

∣

∣

∣

∣

yHFx

yHx

∣

∣

∣

∣

∣

che, data F , e tanto piu grande quanto piu piccolo e |yHx|.Nel caso di autovalore λ di molteplicita algebrica σ(λ) > 1 e molteplicitageometrica τ(λ) si puo arrivare ad una relazione del tipo

|λi(ǫ)− λ| ≤ γ|ǫ|1/µ

per i = 1, . . . , σ(λ), γ costante positiva; se µ > 1 il problema del calcolodi λ puo essere fortemente malcondizionato.

Capitolo 5

Il problema dei MinimiQuadrati

5.1 Le equazioni normali

Sia

Ax = b (5.1)

un sistema lineare con A ∈ IRm×n tale che m ≥ n. Se m > n, il sistema(5.1) ha piu equazioni che incognite e si dice sovradeterminato. Se ilsistema (5.1) non ha soluzione, fissata una norma vettoriale ‖·‖, si cercanoi vettori x ∈ IRn che minimizzano la quantita ‖Ax − b‖. In norma 2, ilproblema diventa quello di determinare un vettore x ∈ IRn tale che

‖Ax− b‖2 = miny∈IRn

‖Ay − b‖2 = µ. (5.2)

Questo problema viene detto problema dei minimi quadrati.

Teorema 5.1 Se A ha rango massimo la soluzione del problema (5.2)esiste, e unica e coincide con la soluzione del sistema lineare

AT Ax = ATb. (5.3)

Tale sistema viene detto sistema delle equazioni normali.

dim. Siano

S(A) = {y ∈ IRm : y = Ax, x ∈ IRn}e

S(A)⊥ = {z ∈ IRm : zTy = 0, ∀y ∈ S(A)}

53

54 CAPITOLO 5. IL PROBLEMA DEI MINIMI QUADRATI

il sottospazio di IRm immagine di A, e il sottospazio ortogonale a S(A). Ilvettore b puo essere cosı decomposto

b = b1 + b2, dove b1 ∈ S(A) e b2 ∈ S(A)⊥

per cui per il residuo

r = b1 − Ax + b2 = y + b2,

dovey = b1 − Ax ∈ S(A) e b2 ∈ S(A)⊥

vale‖r‖22 = (y + b2)

T (y + b2) = ‖y‖22 + ‖b2‖22,in quanto yTb2 = bT

2 y = 0. Poiche solo y dipende da x, si ha che ‖r‖22 eminimo se e solo se b1 = Ax, cioe se e solo se il vettore r ∈ S(A)⊥ ed equindi ortogonale alle colonne di A, cioe

AT r = AT (b− Ax) = 0.

Ne segue quindi che x e soluzione di (5.2) se e solo se e soluzione di (5.3).

5.2 Metodo QR per i minimi quadrati

Sia A di rango massimo e si applichi il metodo di Householder per ottenerela fattorizzazione Q R di A; sara Q ∈ IRm×m ortogonale ed R ∈ IRm×n conla seguente particolare struttura

R =

(

R1

0

)

con R1 ∈ IRn×n triangolare superiore e non singolare essendo H di rangomassimo.Dalla fattorizzazione A = Q R si ha

‖Ax− b‖2 = ‖QRx− b‖2= ‖Q(Rx−QTb)‖2= ‖Rx−QTb‖2= ‖Rx− c‖2

dove c = QTb.Partizionando il vettore c nel modo seguente

c =

(

c1

c2

)

5.2. METODO QR PER I MINIMI QUADRATI 55

con c1 ∈ IRn e c2 ∈ IRm−n, si ha che

Rx− c =

(

R1x− c1

−c2

)

da cuiminx ‖Ax− b‖22 = minx ‖Rx− c‖22

= minx (‖R1x− c1‖22 + ‖c2‖22)= ‖c2‖22 + minx (‖R1x− c1‖22)

Poice R1 e non singolare, esiste una ed una sola soluzione per il sistema

R1x = c1,

sia questa x∗ e sara

minx

(

‖R1x− c1‖22)

= ‖R1x∗ − c1‖22 = 0

Segue che x∗ e la soluzione del problema dei minimi quadrati e

‖c2‖2 = minx‖Ax− b‖2.

56 CAPITOLO 5. IL PROBLEMA DEI MINIMI QUADRATI

Bibliografia

[Bau88] D. Bau, N. Trefethen. Numerical Linear Algebra. SIAM, 1988.

[Bini88] D. Bini, M. Capovani, O. Menchi. Metodi Numerici per l’AlgebraLineare. Zanichelli, 1988.

[Farin88] G. Farin, D. Hansford. The Geometry Toolbox for GraphicsModeling. A.K.Peters, 1988.

57

ARGOMENTI DEL CORSO CALCOLO NUMERICOvvw.web.cs.unibo.it/wiki/images/3/3a/CC-2009_2010-05.pdf1....

Documents

Transcript of ARGOMENTI DEL CORSO CALCOLO NUMERICOvvw.web.cs.unibo.it/wiki/images/3/3a/CC-2009_2010-05.pdf1....