Metodi iterativi per la soluzione di sistemi lineari

Alvise Sommariva

Università degli Studi di Padova

Dipartimento di Matematica

18 aprile 2016

Metodi iterativi e diretti

Problema. Supponiamo che siano A ∈ ℝ^{n×n} (matrice non singolare), b ∈ ℝ^{n} (vettore colonna). Bisogna 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.

Nota. Sia A una matrice invertibile. Allora PA = LU dove

P è una matrice di permutazione,

L è una matrice triangolare inferiore a diagonale unitaria (l_{ii} = 1, ∀i), U una matrice triangolare superiore.

Se Ax = b, nota PA = LU, per risolvere Ax = b da Ax = b ⇔ PAx = Pb ⇔ LUx = Pb si risolve prima Ly = Pb e detta y^{∗} la soluzione, si determina x^{∗} tale che Ux^{∗} = y^{∗}.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  2/ 65

Metodi iterativi e diretti

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. Il costo totale sarà 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, se Ax = b necessariamente Mx − Nx = Ax = b. Di conseguenza da Mx = Nx + b

moltiplicando ambo i membri per M^{−1}, 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.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  4/ 65

Esempio A, E, F

Risulta utile anche scrivere la matrice A come A = D − E − F dove D è una matrice diagonale, E è triangolare inferiore con elementi diagonali nulli, F è triangolare superiore con elementi diagonali nulli.

Definizione Sia A = M − N, con M non singolare. Un metodo iterativo le cui iterazioni sono

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

si dice iterativo stazionario.

La matrice P = M^{−1}N si dice di matrice iterazione.

Nota. Si osservi che la matrice di iterazione M^{−1}N come pure M^{−1}b non dipendono dall’indice di iterazione k.

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

Esempio A, E, F

Esempio

>> 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.
Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  6/ 65

Metodo di Jacobi

Definizione (Metodo di Jacobi, 1846) Sia A = D − E − F ∈ ℝ^{n×n} dove

D è una matrice diagonale, non singolare

E è triangolare inferiore con elementi diagonali nulli, F è triangolare superiore con elementi diagonali nulli.

Sia fissato x^{(0)} ∈ ℝ^{n}. Il metodo di Jacobi genera {x^{(k)}}_{k} con

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

dove

M = D, N = E + F .

Si vede che la matrice di iterazione P = M^{−1}N è in questo caso

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

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  7/ 65

Metodo di Jacobi

Nota. Se D è non singolare il metodo di Jacobi

x^{(k+1)} = ϕ(x^{(k)}) = D^{−1}(E + F )x^{(k)} + D^{−1}b

può essere descritto come metodo delle sostituzioni simultanee

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

Nota. Si osservi che se D è non singolare, ovviamente a_{i,i} ≠ 0 per ogni indice i e quindi la successione è ben definita.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  8/ 65

Metodo di Jacobi

Vediamo perchè questa derivazione. Da Ax = b, con A ∈ ℝ^{n×n} abbiamo

∑_{j<i}a_{i,j}x_{j} + a_{i,i}x_{i} + ∑_{j>i}a_{i,j}x_{j} = ∑_{j}a_{i,j}x_{j} = b_{i}, i = 1, …, n

e quindi evidenziando x_{i}, supposto a_{i,i} ≠ 0 deduciamo

x_{i} = \frac{1}{a_{i,i}} (b_{i} − (∑_{j<i}a_{i,j}x_{j} + ∑_{j>i}a_{i,j}x_{j})), i = 1, …, n

Quindi note delle approssimazioni della soluzione x^{(k)} = {x_{j}^{(k)}} per j = 1, …, n è naturale introdurre il metodo iterativo

x_{i}^{(k+1)} = \frac{1}{a_{i,i}} (b_{i} − (∑_{j<i}a_{i,j}x_{j}^{(k)} + ∑_{j>i}a_{i,j}x_{j}^{(k)})), i = 1, …, n, che si vede facilmente essere il metodo di Jacobi.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  9/ 65

Metodo di Gauss-Seidel

Definizione (Metodo di Gauss-Seidel, 1823-1874) Sia A = D − E − F ∈ ℝ^{n×n} dove

D è una matrice diagonale, non singolare

E è triangolare inferiore con elementi diagonali nulli, F è triangolare superiore con elementi diagonali nulli.

Sia fissato x^{(0)} ∈ ℝ^{n}. Il metodo di Gauss-Seidel genera {x^{(k)}}_{k} con

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

dove

M = D − E, N = F .

Si vede che la matrice di iterazione P = M^{−1}N è in questo caso

P = M^{−1}N = (D − E)^{−1}F .  (3)
Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  10/ 65

Nota. Se D è non singolare il metodo di Gauss-Seidel

x^{(k+1)} = ϕ(x^{(k)}) = (D − E)^{−1}Fx^{(k)} + (D − E)^{−1}b

può essere descritto come metodo delle sostituzioni successive

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

Nota. Si osservi che se D è non singolare, ovviamente a_{i,i} ≠ 0 per ogni indice i e quindi la successione è ben definita.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  11/ 65

Metodo di Gauss-Seidel

Vediamo perchè questa derivazione. Da Ax = b, con A ∈ ℝ^{n×n} abbiamo

∑_{j<i}a_{i,j}x_{j} + a_{i,i}x_{i} + ∑_{j>i}a_{i,j}x_{j} = ∑_{j}a_{i,j}x_{j} = b_{i}, i = 1, …, n

e quindi evidenziando x_{i}, supposto a_{i,i} ≠ 0 deduciamo

x_{i} = \frac{1}{a_{i,i}} (b_{i} − (∑_{j<i}a_{i,j}x_{j} + ∑_{j>i}a_{i,j}x_{j})), i = 1, …, n

Quindi note delle approssimazioni della soluzione x^{(k)} = {x_{j}^{(k)}} per j > i, e già calcolate x^{(k+1)} = {x_{j}^{(k+1)}} per j < i è naturale introdurre il metodo iterativo

x_{i}^{(k+1)} = \frac{1}{a_{i,i}} (b_{i} − (∑_{j<i}a_{i,j}x_{j}^{(k+1)} + ∑_{j>i}a_{i,j}x_{j}^{(k)})), i = 1, …, n.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  12/ 65

Nota. Secondo Forsythe: The Gauss-Seidel method was not known to Gauss and not recommended by Seidel. Secondo Gauss, in una lettera a Gerling (26 dicembre 1823): One could do this even half asleep or one could think of other things during the computations.

SOR

Definizione (Successive over relaxation, (SOR), 1884-1950) Sia A = D − E − F ∈ ℝ^{n×n} dove

D è una matrice diagonale, non singolare

E è triangolare inferiore con elementi diagonali nulli, F è triangolare superiore con elementi diagonali nulli.

Sia fissato x^{(0)} ∈ ℝ^{n} e ω ∈ ℝ\0. Il metodo SOR genera {x^{(k)}}_{k} con

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

dove

M = \frac{D}{ω} − E, N = (\frac{1}{ω} − 1) D + F

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  14/ 65

SOR

Nota. ((cf.[1, p.555])) SOR deriva dallo descrivere esplicitamente una iter. di Gauss-Seidel

z_{i}^{(k+1)} = \frac{1}{a_{ii}} [b_{i} − ∑_{j=1}^{i−1}a_{ij}x_{j}^{(k+1)} − ∑_{j=i+1}^{n}a_{ij}x_{j}^{(k)}]

e sostituirla con la combinazione convessa

x_{i}^{(k+1)} = ωz_{i}^{(k+1)} + (1 − ω)x_{i}^{(k)}.

cioè x_{i}^{(k+1)} = \frac{ω}{a_{ii}} [b_{i} − ∑_{j=1}^{i−1}a_{ij}x_{j}^{(k+1)} − ∑_{j=i+1}^{n}a_{ij}x_{j}^{(k)}] + (1 − ω)x^{(k)}.

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

SOR

Per convincersi dell’equivalenza tra formulazione matriciale ed equazione per equazione, osserviamo che

x_{i}^{(k+1)} = \frac{ω}{a_{ii}} [b_{i} − ∑_{j=1}^{i−1}a_{ij}x_{j}^{(k+1)} − ∑_{j=i+1}^{n}a_{ij}x_{j}^{(k)}] + (1 − ω)x^{(k)} diventa

x^{(k+1)} = ωD^{−1}[b + Ex^{(k+1)} + Fx^{(k)}] + (1 − ω)x^{(k)} e dopo noiosi calcoli in cui si evidenziano x^{(k)} e x^{(k+1)}

\frac{D}{ω} (I − ωD^{−1}E)x^{(k+1)} = \frac{D}{ω} (ωD^{−1}F + (1 − ω))x^{(k)} + b da cui

M = \frac{D}{ω} − E, N = F + \frac{1 − ω}{ω} D = (\frac{1}{ω} − 1) D + F . Nota. Gauss-Seidel è SOR in cui si è scelto ω = 1.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  16/ 65

SOR

Nota. Notiamo che le iterazioni di SOR verificano l’uguaglianza

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

con

M = \frac{D}{ω} − E, N = (\frac{1}{ω} − 1) D + F

ed è

M − N = (\frac{D}{ω} − E) − ((\frac{1}{ω} − 1) D + F ) = −E + D − F = A.

Quindi effettivamente A = M − N.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  17/ 65

I metodi di Richardson stazionari

Definizione (Residuo) Sia A una matrice quadrata n × n, b ∈ ℝ^{n}, x ∈ ℝ^{n}. La quantità

r = b − Ax  (5)

è nota come residuo (relativamente ad A e b). Definizione (Metodo di Richardson stazionario, 1910) Fissato α, un metodo di Richardson stazionario, con matrice di precondizionamento P, verifica

P(x^{(k+1)} − x^{(k)}) = αr^{(k)}.  (6)

dove il residuo alla k-sima iterazione è

r^{(k)} = b − Ax^{(k)}.  (7)
Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  18/ 65

I metodi di Richardson non stazionari

Definizione (Metodo di Richardson non stazionario)

Fissati α_{k} dipendenti dalla iterazione k, un metodo di Richardson non stazionario, con matrice di precondizionamento P, verifica

P(x^{(k+1)} − x^{(k)}) = α_{k}r^{(k)}.  (8)

Nota. Si osservi che se α_{k} = α per ogni k, allora il metodo di Richardson non stazionario diventa stazionario.

I metodi di Richardson

Proposizione.

I metodi di Jacobi e di Gauss-Seidel, SOR, sono metodi di Richardson stazionari.

Dimostrazione. I metodi di Jacobi e di Gauss-Seidel, SOR, sono metodi iterativi del tipo

Mx^{(k+1)} = Nx^{(k)} + b,  (9)

per opportune scelte delle matrici M (che dev’essere invertibile), N tali che A = M − N. Essendo r^{(k)} = b − Ax^{(k)},

M(x^{(k+1)} − x^{(k)}) = Nx^{(k)} + b − Mx^{(k)} = b − Ax^{(k)} = r^{(k)}. (10)

I metodi di Richardson

Ne consegue che i metodi di Jacobi e di Gauss-Seidel, SOR, verificano

M(x^{(k+1)} − x^{(k)}) = r^{(k)}  (11)

In altri termini sono dei metodi di Richardson sono metodi di Richardson stazionari, con α = 1 e matrice di precondizionamento P = M. △

Per quanto riguarda i metodi di Richardson precondizionati e non stazionari, un classico esempio è il metodo del gradiente classico che vedremo in seguito.

Norma di matrici

Definizione (Raggio spettrale) Siano λ_{1}, …, λ_{n} gli autovalori di matrice A ∈ ℝ^{n×n}. La quantità

ρ(A) := max_{i}(|λ_{i}|)

è detto raggio spettrale di A. Definizione (Norma naturale) Sia ∥ · ∥ : ℝ^{n} → ℝ_{+} una norma vettoriale. Definiamo norma naturale (o indotta) di una matrice A ∈ R^{n×n} la quantità

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

Nota. Questa definizione coincide con quella di norma di un operatore lineare e continuo in spazi normati.

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

Esempio ([4], p. 24)

Sia x un arbitrario elemento di R^{n}, A ∈ R^{n×n}.

Norma vettoriale Norma naturale

∥x∥_{1} := ∑_{k=1}^{n} |x_{k}|  ∥A∥_{1} = max_{j} ∑_{i=1}^{n} |a_{i,j}|
∥x∥_{2} := (∑_{k=1}^{n} |x_{k}|^{2})^{1/2} ∥A∥_{2} = ρ^{1/2}(A^{T}A)
∥x∥_{∞} := max_{k} |x_{k}|  ∥A∥_{∞} = max_{i} ∑_{j=1}^{n} |a_{i,j}|

Nota.

Nella norma matriciale 1, di ogni \underline{colonna} si fanno i moduli di ogni componente e li si somma. Quindi si considera il massimo dei valori ottenuti, al variare della riga. Nella norma matriciale ∞, di ogni \underline{riga} si fanno i moduli di ogni componente e li si somma. Quindi si considera il massimo dei valori ottenuti, al variare della colonna.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  23/ 65

Norma di matrici

Per quanto riguarda un esempio chiarificatore in Matlab/Octave

>> A=[1 5 ; 7 1 3 ]

 A =

1  5
7  13
>> norm ( A , 1 )
ans =

18

>> norm ( A , inf )
ans =

20

>> norm ( A , 2 )
ans =

15.5563

>> e i g ( A∗A ’ )
ans =

2 242

>> s q r t (242)
ans =

15.5563

>> raggio_spettrale_A=max( abs ( e i g ( A ) ) )
 raggio_spettrale_A =

15.4261 >>

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

Si dimostra che (cf. [4, p.28]) Teorema Per ogni norma naturale ∥ · ∥ e ogni matrice quadrata A si ha ρ(A) ≤ ∥A∥. Inoltre per ogni matrice A di ordine n e per ogni ϵ > 0 esiste una norma naturale ∥ · ∥ tale che

ρ(A) ≤ ∥A∥ ≤ ρ(A) + ϵ.

e inoltre (cf. [4, p.29], [3, p.232]) Teorema Fissata una norma naturale ∥ · ∥, i seguenti asserti sono equivalenti

1 A^{m} → 0;

2 ∥A^{m}∥ → 0; 3 ρ(A) < 1.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  25/ 65

Sul raggio spettrale

Nota.

Siano {λ_{k}}_{k=1,...,n} autovalori di A ∈ ℝ^{n×n}. Il raggio spettrale

ρ(A) = max_{k}(|λ_{k}|)

non è una norma. Infatti la matrice

( \begin{matrix}0 & 1 \\ 0 & 0 \end{matrix} )

ha raggio spettrale nullo, ma non è la matrice nulla. Osserviamo che dagli esempi il raggio spettrale di una matrice A non coincide in generale con la norma 1, 2, ∞, ma che a volte ρ(A) = ∥A∥_{2} come nel caso di una matrice diagonale A (essendo gli autovalori di una matrice diagonale, proprio i suoi elementi diagonali).

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  26/ 65

Teorema di Hensel

Definizione (Matrice diagonalizzabile) Una matrice A ∈ ℝ^{n×n} è diagonalizzabile se e solo ℝ^{n} possiede una base composta di autovettori di A. Equivalentemente A ∈ ℝ^{n×n} è diagonalizzabile se e solo se è simile ad una matrice diagonale D, cioè esiste S tale che

A = S^{−1}DS.

Definizione (Metodo consistente) Sia A una matrice non singolare e Ax^{∗} = b. Il metodo stazionario

x^{(k+1)} = Px^{(k)} + c

si dice consistente relativamente al problema Ax = b se e solo se

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

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  27/ 65

Lemma Sia x^{(k+1)} = Px^{(k)} + c un metodo iterativo stazionario. Supponiamo sia consistente, relativamente al problema Ax = b, cioè

Ax^{∗} = b implica x^{∗} = Px^{∗} + c.

Allora, posto e^{(k)} = x^{(k)} − x^{(∗)} si ha che e^{(k+m)} = P^{m}e^{(k)}. In particolare, per k = 0 si ottiene e^{(m)} = P^{m}e^{(0)}.

Dimostrazione. Basta osservare che

e^{(k+m)} := x^{(k+m)} − x^{∗} = Px^{(k+m−1)} + c − (Px^{∗} + c)
= P(x^{(k+m−1)} − x^{∗}) = … = P^{m}(x^{(k)} − x^{∗}) = P^{m}e^{(k)}. (12)

Teorema (Teorema di Hensel per matrici diagonalizzabili, 1926)

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

Dimostrazione. Siccome P è diagonalizzabile allora ℝ^{n} possiede una base composta di autovettori di P.

Consideriamo un metodo iterativo stazionario x^{(k+1)} = Px^{(k)} + c in cui scelto x^{(0)} si abbia

e^{(0)} := x^{(0)} − x^{∗} = ∑_{s=1}^{n}c_{s}u_{s}

dove {u_{k}}_{k} è una base di autovettori di P avente autovalori {λ_{k}}_{k}.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  29/ 65

Teorema di Hensel

Supponiamo |λ_{s}| < 1 per s = 1, …, n. Se il metodo è consistente, cioè

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

da un lemma precedente abbiamo e^{(k)} = P^{k}e^{(0)} e da Pu_{s} = λ_{s}u_{s}, x^{(0)} − x^{∗} = ∑_{s=1}^{n} c_{s}u_{s}

P^{k}e^{(0)} = P^{k} (∑_{s=1}^{n}c_{s}u_{s}) = ∑_{s=1}^{n}c_{s}P^{k}u_{s} = ∑_{s=1}^{n}c_{s}P^{k−1}Pu_{s}
= ∑_{s=1}^{n}c_{s}λ_{s}P^{k−1}u_{s} = ∑_{s=1}^{n}c_{s}λ_{s}^{2}P^{k−2}u_{s} = … = ∑_{s=1}^{n}c_{s}λ_{s}^{k}u_{s}.
Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  30/ 65

Dimostrazione. Quindi se |λ_{s}^{k}| < 1 per ogni s = 1, …, n e k = 1, 2, …, abbiamo

∥x^{(k)} − x^{∗}∥ = ∥ ∑_{s=1}^{n}c_{s}λ_{s}^{k}u_{s}∥ ≤ ∑_{s=1}^{n}|c_{s}||λ_{s}^{k}|∥u_{s}∥ → 0.

Se invece per qualche k si ha |λ^{k}| ≥ 1 e c_{k} ≠ 0 allora ∥x^{(k)} − x^{∗}∥ non converge a 0 al crescere di k. Infatti, se λ_{l} ≥ 1 è l’autovalore di massimo modulo, posto

x^{(0)} − x^{∗} = u_{l}

abbiamo da (12)

∥x^{(k)} − x^{∗}∥ = ∥λ_{l}^{k}u_{l}∥ = |λ_{l}^{k}|∥u_{l}∥

che non tende a 0. Di conseguenza non è vero che il metodo è convergente per qualsiasi scelta del vettore x^{(0)}. △

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  31/ 65

Introduciamo il teorema di Hensel nel caso generale.

Teorema (Hensel, 1926, Oldenburger 1950)

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

Nota.

Si noti che il teorema riguarda la convergenza per ogni vettore iniziale x_{0} ed è quindi di convergenza globale. Non si richiede che la matrice P sia diagonalizzabile.

Dimostrazione. ( [3], p. 236) ⇐ Se ρ(P) < 1, allora x = Px + c ha una e una sola sol. x^{∗}. Infatti,

x = Px + c ⇔ (I − P)x = c

e la matrice I − P ha autovalori 1 − λ_{k} con k = 1, …, n tali che

0 < |1 − |λ_{k}|_{ℂ}|_{ℝ} ≤ |1 − λ_{k}|_{ℂ},

poichè |λ_{k}|_{ℂ} ≤ ρ(P) < 1 e quindi

det(I − P) = ∏_{k=1}^{n}(1 − λ_{k}) ≠ 0,

per cui la matrice I − P è invertibile e il sistema (I − P)x = c ha una e una sola soluzione x^{∗}.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  33/ 65

Dimostrazione. Sia e(k) = x^{(k)} − x^{∗}. Come stabilito dal Teorema che lega norme a raggi spettrali, sia inoltre una norma naturale ∥ · ∥ tale che

ρ(P) ≤ ∥P∥ = ρ(P) + (1 − ρ(P))/2 < 1.

Da x^{(k+1)} = Px^{(k)} + c, x = Px + c,

da un lemma precedente e^{(k)} = P^{k}e^{(0)} e quindi essendo ∥ · ∥ una norma naturale

∥e^{(k)}∥ = ∥P^{k}e^{(0)}∥ ≤ ∥P^{k}∥∥e^{(0)}∥.  (13)

Poichè il ρ(P) < 1, abbiamo che ∥P^{k}∥ → 0 da cui per (13) è ∥e^{(k+1)}∥ → 0 e quindi per le proprietà delle norme e^{(k)} → 0 cioè x^{(k)} → x^{∗}.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  34/ 65

Dimostrazione. ⇒ Supponiamo che la successione x^{(k+1)} = Px^{(k)} + c converga a x^{∗} per qualsiasi x^{(0)} ∈ ℝ^{n} ma che sia ρ(P) ≥ 1. Sia λ_{fnsmax} il massimo autovalore in modulo di P e scegliamo x^{(0)} tale che e^{(0)} = x^{(0)} − x^{∗} sia autovettore di P relativamente all’autovalore λ_{fnsmax}. Essendo Pe^{(0)} = λ_{fnsmax}e^{(0)} e e^{(k+1)} = P^{k}e^{(0)} abbiamo che

e^{(k+1)} = λ_{fnsmax}^{k}e^{(0)}

da cui, qualsiasi sia la norma ∥ · ∥, per ogni k = 1, 2, … si ha

∥e^{(k+1)}∥ = |λ_{fnsmax}^{k}|_{ℂ}∥e^{(0)}∥ ≥ ∥e^{(0)}∥

Velocità di convergenza

il che comporta che la successione non è convergente (altrimenti per qualche k sarebbe e^{(k)} < e^{(0)}). △ Proposito.

L’analisi che faremo in questa sezione non è rigorosamente matematica. Ciò nonostante permette di capire il legame tra il raggio spettrale della matrice di iterazione P e la riduzione dell’errore.

Supponiamo di voler approssimare x^{∗} tale che Ax^{∗} = b con il metodo stazionario

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

Supporremo tale metodo consistente e quindi se Ax^{∗} = b allora

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

Inoltre indicheremo con e^{(k)} l’errore

e^{(k)} = x^{(k)} − x^{∗}.
Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  36/ 65

Velocità di convergenza

Da un lemma precedente

e^{(k+m)} = P^{m}e^{(k)}.  (14)

Si dimostra che se A ∈ ℂ^{n×n} e ∥ · ∥ è una norma naturale allora

lim_{k}∥A^{k}∥^{\frac{1}{k}} = ρ(A)

Quindi per k sufficientemente grande si ha ∥P^{k}∥^{1/k} ≈ ρ(P).

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  37/ 65

Sotto queste ipotesi, per k sufficientemente grande, da e^{(k+m)} = P^{m}e^{(k)},

∥P^{k}∥ ≈ ρ^{k}(P),

abbiamo

∥e^{(k+m)}∥ = ∥P^{m}e^{(k)}∥ ≤ ∥P^{m}∥∥e^{(k)}∥ ≈ ρ^{m}(P)∥e^{(k)}∥  (15)

e quindi

∥e^{(k+m)}∥/∥e^{(k)}∥ ⪅ ρ^{m}(P).  (16)

Per cui da

ρ^{m}(P) ≤ ϵ ⇒ ∥e^{(k+m)}∥/∥e^{(k)}∥ ⪅ ρ^{m}(P) ≤ ϵ

applicando il logaritmo naturale ad ambo i membri, per avere una riduzione dell’errore di ϵ, asintoticamente servono al più m iterazioni con,

m log (ρ(P)) ≈ log ϵ

e quindi

m ≈ log ϵ/log (ρ(P)).

Se

R(P) = − log(ρ(P))

è la cosidetta velocità di convergenza asintotica del metodo iterativo relativo a P abbiamo

m ≈ \frac{log ϵ}{log (ρ(P))} = ⌈\frac{− log(ϵ)}{R(P)}⌉ .

Se P è la matrice d’iterazione di un metodo stazionario convergente (e consistente), essendo ρ(P) < 1, minore è ρ(P) necessariamente è maggiore R(P) e si può stimare il numero di iterazioni per ridurre l’errore di un fattore ϵ.

Si desidera quindi cercare metodi metodi con m minore possibile o equivalentemente ρ(P) più piccolo possibile o equivalentemente con R(P) il più grande possibile.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  39/ 65

Alcuni teoremi di convergenza

Definizione (Matrice tridiagonale) 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 & 0 & 9 & 0 \\ 0 & 0 & 4 & 4 & 2 \\ 0 & 0 & 0 & 5 & 3 \end{matrix} )
Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  40/ 65

Alcuni teoremi di convergenza

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

La velocità di convergenza del metodo di Gauss-Seidel è il doppio di quello del metodo di Jacobi.

Nota. Questo significa vuol dire che asintoticamente sono necessarie metà iterazioni del metodo di Gauss-Seidel per ottenere la stessa precisione del metodo di Jacobi.

Alcune definizioni

Definizione (Matrice a predominanza diagonale)

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

La matrice A è a predominanza diagonale per colonne se A^{T} a predominanza diagonale per righe.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  42/ 65

Definizione (Matrice a predom. diagonale in senso stretto)

Se

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

allora la matrice A si dice a predominanza diagonale (per righe) in senso stretto.

La matrice A è a predominanza diagonale (per colonne) in senso stretto se A^{T} a predominanza diagonale per righe 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  43/ 65

Teoremi di convergenza

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

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

Tali teoremi valgono anche se A è a predominanza diagonale per colonne in senso stretto.

Definizione (Matrice simmetrica) La matrice A è simmetrica se A = A^{T}.

Definizione (Matrice definita positiva) Una matrice simmetrica A è definita positiva se ha tutti gli autovalori positivi. Equivalentemente una matrice simmetrica A è definita positiva se

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

Esempio

Ad esempio la matrice

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

è simmetrica e definita positiva. Per convincersene basta vedere che ha tutti gli autovalori positivi.

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

Teorema Sia A simmetrica con elementi diagonali positivi. Allora il metodo di Gauss-Seidel converge se e solo se

0 < w < 2,

A è definita positiva.

Più in generale, Teorema Sia A simmetrica con elementi diagonali positivi. Allora il metodo SOR converge se e solo se

0 < w < 2,

A è definita positiva.

Per una dimostrazione si veda [6, p.215].

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  47/ 65

Test di arresto

Problema. (Test di arresto)

Consideriamo il sistema lineare Ax = b avente un’unica soluzione x^{∗} e supponiamo di risolverlo numericamente con un metodo iterativo stazionario del tipo

x^{(k+1)} = Px^{(k)} + c,

che sia consistente cioè

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

Desideriamo introdurre un test di arresto che interrompa le iterazioni, qualora una certa quantità relativa al sistema lineare Ax = b e alle iterazioni eseguite, sia al di sotto di una tolleranza ϵ > 0 fissata dall’utente.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  48/ 65

Test di arresto: criterio dello step

Proposizione. Sia Ax = b con A non singolare e x^{(k+1)} = Px^{(k)} + c un metodo iterativo stazionario consistente, con ∥P∥_{2} < 1. Allora

∥e^{(k)}∥_{2} := ∥x^{∗} − x^{(k)}∥_{2} ≤ \frac{∥x^{(k+1)} − x^{(k)}∥_{2}}{1 − ∥P∥_{2}}.

Dimostrazione.

Posto Δ^{(k)} := x^{(k+1)} − x^{(k)} e e^{(k)} = x^{∗} − x^{(k)}, essendo e^{(k)} = Pe^{(k−1)} abbiamo

∥e^{(k)}∥_{2} = ∥x^{∗} − x^{(k)}∥_{2} = ∥(x^{∗} − x^{(k+1)}) + (x^{(k+1)} − x^{(k)})∥_{2}
= ∥e^{(k+1)} + Δ^{(k)}∥_{2} = ∥Pe^{(k)} + Δ^{(k)}∥_{2} ≤ ∥P∥_{2} · ∥e^{(k)}∥_{2} + ∥Δ^{(k)}∥_{2}

da cui

(1 − ∥P∥_{2})∥e^{(k)}∥_{2} ≤ ∥Δ^{(k)}∥_{2} ⇔ ∥e^{(k)}∥_{2} := ∥x^{∗} − x^{(k)}∥_{2} ≤ \frac{∥x^{(k+1)} − x^{(k)}∥_{2}}{1 − ∥P∥_{2}}.

Test di arresto: criterio dello step

△ Fissata dall’utente una tolleranza tol, si desidera interrompere il processo iterativo quando

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  49/ 65

∥x^{∗} − x^{(k)}∥ ≤ tol.

Definizione (Criterio dello step)

Il criterio dello step, consiste nell’interrompere il metodo iterativo che genera la successione {x^{(k)}} alla k + 1-sima iterazione qualora

∥x^{(k+1)} − x^{(k)}∥ ≤ tol

(dove tol è una tolleranza fissata dall’utente).

Di seguito desideriamo vedere quando tale criterio risulti attendibile cioè |x^{(k+1)} − x^{(k)}| ≈ |x^{∗} − x^{(k)}|.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  50/ 65

Teorema Sia Ax^{∗} = b con A non singolare e x^{(k+1)} = Px^{(k)} + c un metodo iterativo stazionario

consistente, convergente, con P simmetrica.

Allora

∥e^{(k)}∥_{2} := ∥x^{∗} − x^{(k)}∥_{2} ≤ \frac{∥x^{(k+1)} − x^{(k)}∥_{2}}{1 − ρ(P)}.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  51/ 65

Lemma Se P è simmetrica, allora esistono

una matrice ortogonale U, cioè tale che U^{T} = U^{−1}, una matrice diagonale a coefficienti reali Λ

per cui

P = UΛU^{T}

ed essendo P e Λ simili, le matrici P e Λ hanno gli stessi autovalori {λ_{k}}_{k}.

Dimostrazione. Dal lemma, se P è simmetrica, essendo

∥P∥_{2} := \sqrt{ρ(PP^{T})}, U^{T}U = I

∥P∥_{2} = \sqrt{ρ(PP^{T})} = \sqrt{ρ(UΛU^{T}(UΛU^{T})^{T})}
= \sqrt{ρ(UΛU^{T}UΛU^{T})} = \sqrt{ρ(UΛ^{2}U^{T})}  (17)

Essendo UΛ^{2}U^{T} simile a Λ^{2}, UΛ^{2}U^{T} e Λ^{2} hanno gli stessi autovalori uguali a {λ_{k}^{2}}_{k} e di conseguenza lo stesso raggio spettrale, da cui

ρ(UΛ^{2}U^{T}) = ρ(Λ^{2}).

Da ρ(UΛ^{2}U^{T}) = ρ(Λ^{2}) ricaviamo

∥P∥_{2} = \sqrt{ρ(Λ^{2})} = \sqrt{max_{k}|λ_{k}^{2}|} = \sqrt{max_{k}|λ_{k}|^{2}}
= \sqrt{(max_{k}|λ_{k}|)^{2}} = max_{k}|λ_{k}| = ρ(P).  (18)

Se il metodo stazionario x^{(k+1)} = Px^{(k)} + c è convergente, per il teorema di Hensel ∥P∥_{2} = ρ(P) < 1 e per la precedente proposizione

∥e^{(k)}∥_{2} := ∥x^{∗} − x^{(k)}∥_{2} ≤ \frac{∥x^{(k+1)} − x^{(k)}∥_{2}}{1 − ∥P∥_{2}} = \frac{∥x^{(k+1)} − x^{(k)}∥_{2}}{1 − ρ(P)}.

△

Nota. Il teorema afferma che se la matrice di iterazione P, di un metodo consistente e convergente, è simmetrica, il criterio dello step è affidabile se ρ(P) < 1 è piccolo.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  54/ 65

Test di arresto: criterio del residuo

Definizione (Criterio del residuo)

Il criterio del residuo, consiste nell’interrompere il metodo iterativo che genera la successione {x^{(k)}} alla k + 1-sima iterazione qualora

∥r^{(k)}∥ ≤ tol

dove tol è una tolleranza fissata dall’utente, r^{(k)} = b − Ax^{(k)}.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  55/ 65

Test di arresto: criterio del residuo

Teorema Sia Ax^{∗} = b con A non singolare. Sia {x^{(k)}} la successione generata da un metodo iterativo. Allora

\frac{∥x^{(k)} − x^{∗}∥}{∥x^{∗}∥} = \frac{∥e^{(k)}∥}{∥x^{∗}∥} ≤ κ(A)\frac{∥r^{(k)}∥}{∥b∥}.

dove κ(A) è il numero di condizionamento di A.

Dimostrazione. Se b = Ax^{∗} allora r^{(k)} = Ae^{(k)} in quanto

r^{(k)}b − Ax^{(k)} = Ax^{∗} − Ax^{(k)} = A(x^{∗} − x^{(k)}) = Ae^{(k)}.

Osserviamo che

essendo A invertibile e r^{(k)} = Ae^{(k)}, si ha e^{(k)} = A^{−1}r^{(k)} e

∥e^{(k)}∥ = ∥A^{−1}r^{(k)}∥ ≤ ∥A^{−1}∥∥r^{(k)}∥;

poichè b = Ax^{∗} abbiamo ∥b∥ ≤ ∥A∥∥x^{∗}∥ e quindi \frac{1}{∥x^{∗}∥} ≤ \frac{∥A∥}{∥b∥}.

Test di arresto: criterio del residuo

Di conseguenza, denotato con κ(A) = ∥A∥∥A^{−1}∥ il numero di condizionamento (necessariamente maggiore o uguale a 1), da

∥e^{(k)}∥ = ∥A^{−1}r^{(k)}∥ ≤ ∥A^{−1}∥∥r^{(k)}∥

\frac{1}{∥x^{∗}∥} ≤ \frac{∥A∥}{∥b∥}

se x^{∗} ≠ 0 abbiamo

\frac{∥e^{(k)}∥}{∥x^{∗}∥} ≤ \frac{∥A∥}{∥b∥}∥e^{(k)}∥ ≤ \frac{∥A∥}{∥b∥} · ∥A^{−1}∥∥r^{(k)}∥ ≤ κ(A) \frac{∥r^{(k)}∥}{∥b∥}

Quindi

\frac{∥x^{(k)} − x^{∗}∥}{∥x^{∗}∥} = \frac{∥e^{(k)}∥}{∥x^{∗}∥} ≤ κ(A) \frac{∥r^{(k)}∥}{∥b∥} ≤ tol.

△

Nota. Dal teorema si evince che il criterio \frac{∥r^{(k)}∥}{∥b∥} ≤ tol può essere conservativo per κ(A) ≫ 1.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  57/ 65

I 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 dell’energia

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

in quanto

grad(ϕ(x)) = Ax − b = 0 ⇔ Ax = b.

Quindi invece di calcolare la soluzione del sistema lineare, intendiamo calcolare il minimo del funzionale ϕ.

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.

Definizione (Metodo del gradiente classico, Cauchy 1847)

Il metodo del gradiente classico è un metodo di discesa x^{(k+1)} = x^{(k)} + α_{k}p^{(k)}

con α_{k} e p^{(k)} scelti così da ottenere la massima riduzione del funzionale dell’energia a partire dal punto x^{(k)}. Differenziando ϕ(x) = \frac{1}{2}x^{T}Ax − b^{T}x, x ∈ ℝ^{n} si mostra che tale scelta coincide con lo scegliere

p^{(k)} = r^{(k)}, α_{k} = \frac{∥r^{(k)}∥_{2}^{2}}{(r^{(k)})^{T}Ar^{(k)}}.  (19)

Nota. Con qualche facile conto si vede che è un metodo di Richardson non stazionario con P = I e parametro α_{k} definito da (19).

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  60/ 65

Il metodo del gradiente coniugato.

Definizione (Gradiente coniugato, Hestenes-Stiefel, 1952) Il metodo del gradiente coniugato è un metodo di discesa

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

con

α_{k} = \frac{(r^{(k)})^{T}r^{(k)}}{(p^{(k)})^{T}Ap^{(k)}},
p^{(0)} = r^{(0)},
p^{(k)} = r^{(k)} + β_{k}p^{(k−1)}, k = 1, 2, …
β_{k} = \frac{(r^{(k)})^{T}r^{(k)}}{(r^{(k−1)})^{T}r^{(k−1)}} = \frac{∥r^{(k)}∥_{2}^{2}}{∥r^{(k−1)}∥_{2}^{2}}.

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  61/ 65

Il metodo del gradiente coniugato: proprietà

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

Teorema Se A è una matrice simmetrica e definita positiva di ordine n, allora il metodo del gradiente coniugato è convergente e fornisce in aritmetica esatta la soluzione del sistema Ax = b in al massimo n iterazioni.

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

Il metodo del gradiente coniugato: proprietà

Definizione (Spazio di Krylov)

Lo spazio

K_{k} = span(r^{(0)}, Ar^{(0)}, …, A^{k−1}r^{(0)})

per k ≥ 1 si dice spazio di Krylov.

Teorema Sia

K_{k} = span(r^{(0)}, Ar^{(0)}, …, A^{k−1}r^{(0)})

per k ≥ 1. Allora la k-sima iterata dal metodo del gradiente coniugato, minimizza il funzionale ϕ nell’insieme x^{(0)} + K_{k}.

Per una dimostrazione si veda [7, p.12]. Si osservi che essendo la k-sima iterazione del gradiente classico pure in x^{(0)} + K_{k}, il gradiente classico non minimizza in generale il funzionale ϕ in x^{(0)} + K_{k}.

Alvise Sommariva  Metodi iterativi per la soluzione di sistemi lineari  63/ 65

Il metodo del gradiente coniugato: proprietà

Teorema Se

A è simmetrica e definita positiva, ∥x∥_{A} = \sqrt{x^{T}Ax}, e_{k} = x^{∗} − x^{(k)},

allora

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

Nota. 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}|}, \{λ_{i}\} = spettro(A).

Bibliografia

K. Atkinson, Introduction to Numerical Analysis, Wiley, 1989. K. Atkinson e W. Han, Theoretical Numerical Analysis, Springer, 2001. D. Bini, M. Capovani e O. Menchi, Metodi numerici per l’algebra lineare, Zanichelli, 1988. V. Comincioli, Analisi Numerica, metodi modelli applicazioni, Mc Graw-Hill, 1990. S.D. Conte e C. de Boor, Elementary Numerical Analysis, 3rd Edition, Mc Graw-Hill, 1980. L.A. Hageman e D.M. Young, Applied Iterative Methods, Dover, 2004. C.T. Kelley, Iterative Methods for Linear and Nonlinear Equations, SIAM, 1995. The MathWorks Inc., Numerical Computing with Matlab, http://www.mathworks.com/moler. Netlib, http://www.netlib.org/templates/matlab/. A. Quarteroni e F. Saleri, Introduzione al calcolo scientifico, Springer Verlag, 2006.