Metodi iterativi per la soluzione di sistemi lineari

Alvise Sommariva

Università degli Studi di Padova Dipartimento di Matematica Pura e Applicata

12 gennaio 2012

Metodi iterativi e diretti

Supponiamo che siano A ∈ ℝ^{n×n} (matrice non singolare), b ∈ ℝ^{n} (vettore colonna) e di dover risolvere il problema Ax = b avente sol. unica x^{∗}. A tal proposito si può utilizzare la fatt. LU con pivoting. Il costo computazionale è in generale di O(n^{3}/3) operazioni moltiplicative. Questo diventa proibitivo se n è particolarmente elevato.

L’idea dei metodi iterativi è quello di ottenere una successione di vettori x^{(k)} → x^{∗} cosicchè per k̄ ≪ n sia x^{(k̄)} ≈ x^{∗}. In generale, la soluzione non è ottenuta esattamente come nei metodi diretti in un numero finito di operazioni (in aritmetica esatta), ma quale limite, accontentandosi di poche iterazioni ognuna dal costo quadratico. Quindi il costo totale sarà di ordine O(k̄ · n^{2}).

Metodi iterativi stazionari

Supponiamo che la matrice non singolare A ∈ ℝ^{n×n} sia tale che

A = M − N, con M non singolare.

Allora Mx − Nx = Ax = b e quindi Mx = Nx + b. Moltiplicando ambo i membri per M^{−1} e posto ϕ(x) = M^{−1}Nx + b abbiamo x = M^{−1}Nx + M^{−1}b = ϕ(x). Viene quindi naturale utilizzare la succ. del metodo di punto fisso

x^{(k+1)} = ϕ(x^{(k)}) = M^{−1}Nx^{(k)} + M^{−1}b

La matrice P = M^{−1}N si dice di iterazione e non dipende, come pure b dall’indice di iterazione k. Per questo motivo tali metodi si chiamano iterativi stazionari.

Quale utile notazione, sia inoltre A = D − E − F con D la matrice diagonale estratta da A, E, F rispettivemente triangolari inferiore e superiore.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  3/ 24

Metodo di Jacobi

Il metodo di Jacobi fu scoperto nel 1845, nell’ambito di alcune ricerche su problemi di piccole oscillazioni che comportavano alla risoluzione di sistemi lineari con matrici diagonalmente dominanti. Nel caso del metodo di Jacobi si ha

M = D, N = E + F  (1)

e quindi

P = M^{−1}N = D^{−1}(E + F ) = D^{−1}(D − D + E + F )
= D^{−1}(D − A) = I − D^{−1}A  (2)

Si osservi che se D è non singolare allora il metodo di Jacobi, almeno in questa versione di base, non può essere utilizzato visto che in (5) non ha senso la scrittura D^{−1}.

Metodo di Jacobi

Qualora sia a_{ii} ≠ 0 per ogni i = 1, …, n, il metodo di Jacobi può essere descritto come

x_{i}^{(k+1)} = (b_{i} − ∑_{j=1}^{i−1}a_{ij}x_{j}^{(k)} − ∑_{j=i+1}^{n}a_{ij}x_{j}^{(k)})/a_{ii}, i = 1, …, n. (3)

L’idea è che se Ax = B allora ∑_{j=1}^{n} a_{s,j}x_{j} = b allora

a_{i,i}x_{i} = b − ∑_{j=1}^{i−1}a_{s,j}x_{j} − ∑_{j=i+1}^{n}a_{s,j}x_{j}

da cui (3) ricordando il metodo delle appross. successive, qualora a_{i,i} ≠ 0. Si osservi che ogni componente di x^{(k+1)} dipende dai soli valori all’iterazione precedente x^{(k)}, e questo rende il metodo facilmente implementabile in computers paralleli. Da (3) si capisce perchè tale metodo è noto anche come metodo delle sostituzioni simultanee.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  5/ 24

Metodo di Gauss-Seidel

Il metodo di Gauss-Seidel fu scoperto nel 1874, da studi preliminari di Gauss (1823) completati dal suo allievo Seidel per lo studio di problemi ai minimi quadrati del tipo Sx = f con S non quadrata, che venivano risolti quali soluzione del sistema di equazioni normali S^{T}Sx = S^{T}f . Mentre Gauss oltre a problemi di Astronomia era interessato a problemi di Geodesia (triangolazione di Hannover usando una catena di 26 triangoli), Seidel si interessava alla risoluzione di un sistema di equazioni con 72 incognite per uno studio di luminosità stellare. Il metodo di Gauss-Seidel è definito quale metodo stazionario in cui

M = D − E, N = F  (4)

e quindi

P = M^{−1}N = (D − E)^{−1}F  (5)

Metodo di Gauss-Seidel

Similmente al metodo di Jacobi, possiamo riscrivere più semplicemente anche Gauss-Seidel come

x_{i}^{(k+1)} = (b_{i} − ∑_{j=1}^{i−1}a_{ij}x_{j}^{(k+1)} − ∑_{j=i+1}^{n}a_{ij}x_{j}^{(k)})/a_{ii}.  (6)

Da (6) si capisce perchè tale metodo è noto anche come metodo delle sostituzioni successive.

Alcune definizioni

▶ Siano {λ_{k}}_{k=1,...,n} gli autovalori di una matrice A. Si definisce raggio spettrale ρ(A) la quantità

max_{k=1,...,n}|λ_{k}|.

▶ Sia ∥ · ∥ : ℝ^{n} → ℝ_{+} una norma vettoriale. Definiamo norma naturale di una matrice A ∈ R^{n×n} la quantità

∥A∥ := sup_{x∈ℝ^{n},x≠0} \frac{∥Ax∥}{∥x∥}.

Per ogni norma naturale ∥ · ∥ si ha ρ(A) ≤ ∥A∥. ▶ Un metodo iterativo stazionario x^{(k+1)} = Px^{(k)} + c è

consistente per la soluzione di Ax = b se per x^{∗} tale che Ax^{∗} = b si ha x^{∗} = Px^{∗} + c.

Una matrice A si dice diagonalizzabile se vale uno qualsiasi di queste (equivalenti) affermazioni:

▶ è simile a una matrice diagonale, cioè A = S^{−1}DS con D

matrice diagonale. ▶ esiste una base di autovettori {u_{k}} (con autovalori {λ_{k}}). ▶ tutti gli autovalori di A sono distinti;

NB: Una matrice A simmetrica (cioè A = A^{T}) è diagonalizzabile.

Teorema di Hensel

Teorema. Se P è diagonalizz. allora un metodo iter. staz. cons. x^{(k+1)} = Px^{(k)} + c converge per ogni vett. iniziale x_{0} ↔ ρ(P) < 1.

(⇐) Se la matrice P è diagonalizz. allora esiste una base di autovett. {u_{k}} con autovalori {λ_{k}}. Poichè il metodo è consistente

x^{∗} = Px^{∗} + c.

Scelto arbitrariamente x^{0}, abbiamo per certi c_{1}, …, c_{n} che

x^{0} − x^{∗} = ∑_{s=1}^{n}c_{s}u_{s}

e quindi da x^{(k)} = Px^{(k−1)} + c e x^{∗} = Px^{∗} + c, per linearità

x^{(k)}−x^{∗} = Px^{(k−1)}+c−(Px^{∗}+c) = Px^{(k−1)}−Px^{∗} = P(x^{(k−1)}−x^{∗}).

Teorema di Hensel

Ma allora, ripetendo più volte tale ragionamento

x^{(k)} − x^{∗} = P(x^{(k−1)} − x^{∗}) = P(P(x^{(k−2)} − x^{∗}))
= P^{2}(x^{(k−2)} − x^{∗}) = … = P^{k}(x^{0} − x^{∗}).

Quindi x^{(k)} − x^{∗} = P^{k}(x^{0} − x^{∗}) = P^{k}(∑_{s=1}^{n}c_{s}u_{s}) = ∑_{s=1}^{n}c_{s}P^{k}u_{s} = ∑_{s=1}^{n}c_{s}λ_{s}^{k}u_{s}

Così da |λ_{s}^{k}| ≤ (ρ(P))^{k} < 1 per ogni s = 1, …, n e k = 1, 2, …,

∥x^{(k)} − x^{∗}∥ = ∥ ∑_{s=1}^{n}c_{s}λ_{s}^{k}u_{s}∥ ≤ ∑_{s=1}^{n}|c_{s}||λ_{s}^{k}|∥u_{s}∥
≤ (ρ(P))^{k} ∑_{s=1}^{n}|c_{s}|∥u_{s}∥ → 0
Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  11/ 24

(⇒) Supponiamo per assurdo che per qualche k sia |λ_{s}| ≥ 1. Posto c ≠ 0 sia x^{(0)} = x^{∗} + cu_{s} con Pu_{s} = λ_{s}u_{s}. Similmente a quanto visto in precedenza:

∥x^{(k)} − x^{∗}∥ = ∥cλ_{s}^{k}u_{s}∥ ≤ |λ_{s}^{k}|∥cu_{s}∥ → ∞

Di conseguenza abbiamo trovato un punto iniziale x^{(0)} per cui il metodo non converge, il che è assurdo.

Esiste una versione più generale del precedente teorema (per P non diagonalizzabili).

Teorema. Un metodo iterativo stazionario consistente x^{(k+1)} = Px^{(k)} + c converge per ogni vettore iniziale x_{0} se e solo se ρ(P) < 1.

Alcuni teoremi di convergenza

Una matrice A si dice tridiagonale se A_{i,j} = 0 qualora |i − j| > 1.

Esempio:

( \begin{matrix}1 & 2 & 0 & 0 & 0 \\ 1 & 5 & 7 & 0 & 0 \\ 0 & 2 & 6 & 9 & 0 \\ 0 & 0 & 4 & 4 & 2 \\ 0 & 0 & 0 & 5 & 3 \end{matrix} )

Teorema. Per matrici tridiagonali A = (a_{i,j}) con componenti diagonali non nulle, i metodi di Jacobi e Gauss-Seidel sono o entrambi convergenti o divergenti e il tasso di convergenza del metodo di Gauss-Seidel è il doppio di quello del metodo di Jacobi (il che vuol dire che asintoticamente sono necessarie metà iterazioni del metodo di Gauss-Seidel per ottenere la stessa precisione del metodo di Jacobi).

Ricordiamo che A è a predominanza diagonale (per righe) se per ogni i = 1, …, n risulta

|a_{i,i}| ≥ ∑_{j=1,j≠s}^{n} |a_{i,j}|

e per almeno un indice s si abbia

|a_{s,s}| > ∑_{j=1,j≠s}^{n} |a_{s,j}|.

Se per tutti gli indici s si ha |a_{s,s}| > ∑_{j=1,j≠s}^{n} |a_{s,j}| allora la predominanza diagonale è in senso stretto. Ad esempio la matrice

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

è a predominanza diagonale per righe (non stretta).

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  15/ 24

Teoremi di convergenza

▶ Sia A una matrice quadrata a predominanza diagonale per righe, in senso stretto. Allora il metodo di Jacobi converge alla soluzione di Ax = b, qualsiasi sia il punto x^{(0)} iniziale.

▶ Sia A è a predominanza diagonale per righe, in senso stretto. Allora il metodo di Gauss-Seidel risulta convergente, qualsiasi sia il punto x^{(0)} iniziale.

Tali teoremi valgono se A^{T} è a predominanza diagonale per righe in senso stretto (cioè A è a predominanza diagonale per colonne in senso stretto).

La matrice A è simmetrica se A = A^{T}. Una matrice A è definita positiva se ha tutti gli autovalori positivi. Ad esempio la matrice

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

è simmetrica e definita positiva:

>> A=[4 −1 0 ; −1 4 −1; 0 −1 4 ] ;
>> e i g ( A ) % AUTOVALORI DI A .
ans =

2.5858 4.0000 5.4142

>>

Alcuni teoremi di convergenza

Ricordiamo che equivalentemente una matrice A è definita positiva se

x^{T}Ax > 0, per ogni x ≠ 0.

Teorema. Sia A una matrice simmetrica, non singolare con elementi principali a_{i,i} ≠ 0. Allora il metodo di Gauss-Seidel è convergente per qualsiasi scelta del punto iniziale x^{(0)} se e solo se A è definita positiva.

I metodi di discesa

Una classica famiglia di metodi iterativi è quella dei metodi di discesa. Sia A una matrice simmetrica definita positiva. Si osserva che se x^{∗} è l’unica soluzione di Ax = b allora è pure il minimo del funzionale

ϕ(x) = \frac{1}{2}x^{T}Ax − b^{T}x, x ∈ ℝ^{n}

Un generico metodo di discesa consiste nel generare una successione

x^{(k+1)} = x^{(k)} + α_{k}p^{(k)}

dove p^{(k)} è una direzione fissata secondo qualche criterio. Lo scopo ovviamente è che

ϕ(x^{(k+1)}) < ϕ(x^{(k)}),

e che il punto x^{∗}, in cui si ha il minimo di ϕ, venga calcolato rapidamente. Il gradiente coniugato è un classico esempio di metodo di discesa.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  19/ 24

Il metodo del gradiente coniugato.

Il metodo del grad. coniugato fu descritto nel 1952 da Hestenes e Stiefel ma non venne molto utilizzato fino al 1971, quando Reid suggerì il suo utilizzo per la risoluzione di sistemi sparsi (cioè con molte componenti nulle) di grandi dimensioni. Il metodo non è iterativo stazionario. Supponiamo di dover risolvere il sistema lineare Ax = b. Con il termine residuo in x^{(k)} si intende la quantità r^{(k)} = b − Ax^{(k)}. La succ. delle iterazioni del gradiente coniugato è quella propria dei metodi di discesa,

x^{(k+1)} = x^{(k)} + α_{k}p^{(k)}, α_{k} = \frac{(r^{(k)})^{T}r^{(k)}}{(p^{(k)})^{T}Ap^{(k)}}

dove p^{(0)} = r^{(0)} e

p^{(k)} = r^{(k)} + β_{k}p^{(k−1)}, β_{k} = \frac{(r^{(k)})^{T}r^{(k)}}{(r^{(k−1)})^{T}r^{(k−1)}}.

Con questa scelta si prova che p^{(k)} e p^{(k−1)} sono A-coniugati.

(p^{(k)})^{T}Ap^{(k−1)} = 0.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  20/ 24

Il metodo del gradiente coniugato: proprietà

Il metodo del gradiente coniugato ha molte proprietà particolari. Ne citiamo alcune.

▶ Se A è una matrice simmetrica e definita positiva di ordine n, si può dimostrare che il metodo è convergente e fornisce in aritmetica esatta la soluzione del sistema Ax = b in al massimo n iterazioni.

Questo teorema tradisce un po’ le attese, sia perchè in generale i calcoli non sono compiuti in aritmetica esatta, sia perchè in molti casi della modellistica matematica n risulta essere molto alto.

Il metodo del gradiente coniugato: proprietà

Si può dimostrare che se A è simmetrica e definita positiva,

∥x∥_{A} = \sqrt{x^{T}Ax}

e

e_{k} = x^{∗} − x^{(k)}

allora

∥e_{k}∥_{A} ≤ (\frac{\sqrt{K_{2}(A)} − 1}{\sqrt{K_{2}(A)} + 1})^{2k} ∥e_{0}∥_{A}.

Questo risultato stabilisce che la convergenza del gradiente coniugato è lenta qualora si abbiano alti numeri di condizionamento

K_{2}(A) := ∥A∥_{2}∥A^{−1}∥_{2} = \frac{max_{i} |λ_{i}|}{min_{j} |λ_{j}|}

(ove al solito {λ_{i}} sono gli autovalori di A).

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  22/ 24

Gradiente coniugato in Matlab

Per quanto riguarda il codice del Gradiente Coniugato, un esempio è il file cg.m tratto da Netlib:

f u n c t i o n [ x , e r r o r , iter , f l a g ]=cg ( A , x , b , M , max_it , tol )
f l a g = 0 ; iter = 0 ; bnrm2 = norm ( b ) ;
i f ( bnrm2 = 0 . 0 ) , bnrm2 = 1 . 0 ; end
r = b − A∗x ; e r r o r = norm ( r ) / bnrm2 ;
i f ( e r r o r < tol ) r e t u r n , end

f o r iter = 1 : max_it

z = M \ r ; rho = ( r ’ ∗ z ) ;
i f ( iter > 1 )
beta = rho / rho_1 ; p = z + beta ∗p ;

e l s e

    p = z ;

end

  q = A∗p ; alpha = rho / ( p ’ ∗ q ) ; x = x + alpha ∗ p ;

r = r − alpha∗q ; e r r o r = norm ( r ) / bnrm2 ;
i f ( e r r o r <= tol ) , break , end

  rho_1 = rho ;

end

i f ( e r r o r > tol ) f l a g = 1 ; end
Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  23/ 24

Esercizio in Matlab

▶ Ha Matlab una sua versione del gradiente coniugato? Se la risposta è negativa cercare in internet una versione del gradiente coniugato per Matlab.

▶ Matlab dispone di una gallery di matrici. Aiutandosi con

>> h e l p g a l l e r y

si calcoli la matrice di Poisson P_{20} di ordine 20. ▶ Sia b il vettore composto di componenti uguali a 1, avente lo

stesso numero di righe di P_{20}. Si risolva col gradiente coniugato il problema P_{20}x = b. Quante iterazioni servono?