Metodi iterativi per la soluzione di sistemi lineari

Alvise Sommariva

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

16 maggio 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/ 17

Esempio A, E, F

>> A=[1 2 3 4 ; 5 6 7 2 ; 8 9 1 2 ; 3 4 5 1 ]

A =

1  2  3  4
5  6  7  2
8  9  1  2
3  4  5  1
>> E=−( t r i l ( A )−d i a g ( d i a g ( A ) ) )

E =

0  0  0  0
−5  0  0  0
−8  −9  0  0
−3  −4  −5  0
>> F=−( t r i u ( A )−d i a g ( d i a g ( A ) ) )

F =

0  −2  −3  −4
0  0  −7  −2
0  0  0  −2
0  0  0  0
>> % A=d i a g ( d i a g (A) )−E−F .

Metodo di Jacobi e di Gauss-Seidel

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 (4) non ha senso la scrittura D^{−1}. Il metodo di Gauss-Seidel è definito quale metodo stazionario in cui

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

e quindi

P = M^{−1}N = (D − E)^{−1}F  (4)
Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  5/ 17

Teorema di Hensel

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

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

Teorema. Un metodo iterativo stazionario consistente x^{(k+1)} = Px^{(k)} + c converge per ogni vettore iniziale x_{0} e termine noto c 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).

Alcune definizioni

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  8/ 17

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

Alcune definizioni

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  12/ 17

Il metodo del gradiente coniugato (1952).

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.

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  15/ 17

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  16/ 17

Esercizio in Matlab I

▶ Eseguire una routine iterstat che implementi, date due matrici M ed N, un metodo iterativo stazionario. Quale criterio di arresto si utilizzi

∥x^{(k+1)} − x^{(k)}∥_{∞} < tol

con tol tolleranza richiesta dall’utente (es. tol = 10^{(−12)}). ▶ Usare iterstat per definire jacobi e gs che implementano

il metodo di Jacobi e Gauss Seidel, per risolvere il problema Ax = b. Se necessario usare i comandi diag, tril e triu.

▶ Si calcoli la matrice di Poisson P_{20} di ordine 20 aiutandosi con

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

Che dimensione ha la matrice P_{20}?

Esercizio in Matlab II

▶ 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? E quante con Jacobi e Gauss-Seidel? Come si puó controllare che e’ simmetrica attraverso la norma di A − A^{T}? E vedere che é definita positiva? E che il raggio spettrale é minore di 1?