I METODI DI DISCESA

I METODI DI DISCESA
1. Generalità sui metodi di discesa
Per la risoluzione di un sistema lineare Ax = b con matrice A reale, simmetrica e
definita positiva, un’altra famiglia di metodi iterativi è data dai così detti metodi di
discesa.
Il risultato teorico alla base di questi metodi è il seguente:
Teorema:
Sia A ∈ Rnxn matrice simmetrica e definita positiva, allora la soluzione del sistema
Ax = b
lineare
(1)
coincide con il punto di minimo della seguente funzione quadratica
n
1
1 T
1 n n
T
F ( x ) = ( Ax , x) − (b, x) = x Ax − b x = ∑∑ aij x i x j − ∑ bi xi
2
2
2 i =1 j =1
i =1
(2)
ove la forma quadratica
Q(x)=(Ax,x)=xTAx è positiva per x≠ 0.
Il teorema afferma quindi che risolvere il sistema (1) equivale a minimizzare la
funzione (2).
Dimostrazione:
Se consideriamo il sistema (1) e sostituiamo in esso al posto di x un vettore tentativo
x(k) si otterrà un vettore residuo
r(k) = Ax(k) - b
di componenti r1( k ) , r2( k ) ,...,rn( k ) .
La soluzione del sistema si otterrà modificando successivamente il vettore x(k) in
modo tale che il vettore r(k) diventi sempre più piccolo fino a scomparire; si ha cioè
che il vettore x* corrisponde a r = Ax* - b = 0.
1
Ora consideriamo la funzione quadratica:
1
F ( x ) = ( Ax , x) − (b , x )
2
e cerchiamo il suo punto di minimo.
A questo scopo deriviamo parzialmente la funzione quadratica ed eguagliamo a zero
le derivate; per il fatto che la Q(x) è positiva il punto trovato sarà proprio il minimo
della F.
∂F n
= ∑ aij x j − bi = 0 i = 1,.., n
∂xi j =1
ma questo equivale a scrivere
r = ∇ F = Ax*- b = 0.
Quindi il vettore x* che minimizza la funzione F(x) coincide con la soluzione del
sistema lineare (1); viceversa la soluzione del sistema (1) con matrice A simmetrica
definita positiva minimizza la corrispondente funzione quadratica (2).
Il risultato di questo teorema ci permette di affermare che per la risoluzione di sistemi
lineari con matrice simmetrica definita positiva in generale possono essere usati i
metodi per determinare il minimo di una funzione quadratica, noti come metodi di
discesa.
Questi metodi consistono nel determinare, a partire da un vettore iniziale x, una
direzione p opportuna e nel correggere x in questa direzione in modo che il valore
della funzione quadratica F nel nuovo punto x + t⋅p diminuisca, cioè
F ( x + tp) < F ( x) .
Perché ciò avvenga il parametro t e la direzione p devono essere scelti in modo
opportuno. I differenti metodi di discesa sono caratterizzati dalla scelta della
direzione di discesa p fra le direzioni di discesa ammissibili.
2
La determinazione del parametro t
che permette di rendere minima la F nella
direzione p è quello che differenzia i metodi suddetti quando vengono impiegati per
la soluzione del sistema lineare Ax = b (quindi in cui A sia nota) da quando vengono
impiegati per determinare il minimo di una funzione qualsiasi (in cui quindi la
matrice A è la matrice delle derivate parziali seconde e quindi non è nota).
2. Scelta del parametro t
Nel caso in cui la matrice A sia nota si vede infatti che
1
F ( x + tp ) = ( A( x + tp), ( x + tp)) − (b, ( x + tp ))
2
1
1
= ( Ax , x ) + t ( Ax , p ) + t 2 ( Ap , p ) − (b, x ) − t (b, p )
2
2
1
= F ( x ) + t 2 ( Ap , p ) + t (r , p ),
r = Ax − b,
2
è una funzione quadratica in t.
Per determinarne il minimo rispetto a t basta considerare
dF
= t ( Ap , p ) + ( r, p )
dt
e uguagliarlo a 0; si ottiene per il parametro t il valore
tmin
( r , p)
rT p
=−
= T
( Ap, p ) p Ap
(3)
Questo significa che, una volta determinata la direzione p, ovvero scelto il metodo di
minimizzazione, la relazione (3) ci fornisce il valore da assegnare al parametro t per
ottenere il minimo valore possibile della F lungo la direzione p.
3
Osservazione: il valore di t fornito dalla (3) ci fornisce un minimo per F in quanto se
effettuiamo la derivata seconda in t otteniamo
d 2F
= ( Ap , p)
2
dt
Che è una quantità positiva per ogni direzione di p.
Nel caso in cui il metodo venga applicato per determinare il minimo di una funzione
e quindi la matrice A non sia nota, il valore di t ottimale si ottiene risolvendo un
problema di minimo monodimensionale per la funzione F(t).
Nel punto di minimo ottenuto usando la formula (3) per il valore del parametro t
abbiamo il seguente risultato:
Teorema:
Nel punto di minimo x’ = x + tmin p, ottenuto muovendosi lungo la direzione p con t
dato dalla (3), il vettore residuo r’ = Ax’ – b risulta ortogonale alla direzione p, cioè
(r ', p ) = 0
(*)
Dimostrazione:
Essendo
r ' = Ax '−b = A ( x + tp ) − b = r + tAp
Si ottiene
(r ' , p) = ( r, p) + t ( Ap, p)
che risulta 0 scegliendo t secondo la (3).
4
3. Condizioni di ammissibilità per la direzione di discesa
Per quanto riguarda la scelta delle direzioni di discesa p se si considera la relazione
t min = −
(r , p )
( Ap , p )
si può affermare che p non deve essere ortogonale al residuo r, o, in modo
equivalente, non deve essere ortogonale al gradiente di F (∇F) perché questo
porterebbe a tmin=0.
Se inoltre consideriamo lo sviluppo in serie di Taylor della F, cioè
F ( x + tp ) = F ( x ) + t ⋅ (∇F , p ) + ...
e richiediamo che nel punto (x + tp) si abbia
F ( x + tp ) < F ( x )
per t ≥ 0 ,
(4)
allora la direzione p deve soddisfare la seguente condizione
(∇F , p ) < 0
che rappresenta la condizione di direzione ammissibile.
Questa condizione ci dice che l’angolo fra la nuova direzione di discesa e il gradiente
della F(x) deve avere coseno negativo (ricordiamo che (∇F , p ) =| ∇F | ⋅ | p | cosθ
dove l’angolo θ è l’angolo fra i due vettori) cioè l’angolo θ deve essere maggiore di
π/2.
Poiché la condizione di ammissibilità per la direzione di discesa è data in funzione
del gradiente della F(x) i metodi di discesa vengono anche chiamati metodi del
gradiente.
Es: Interpretazione geometrica dei metodi di discesa
Nel caso n = 2 la funzione F(x)=cost è rappresentata da ellissi concentriche il cui
centro coincide con il minimo della funzione quadratica F(x) e costituisce la
5
soluzione del problema. Il seguente grafico mostra quanto affermato dalla condizione
di ammissibilità per la direzione di discesa.
4. Metodo della Discesa più Ripida (Steepest Descent)
Il metodo di Discesa più Ripida (Steepest Descent) è caratterizzato dalla scelta,
ad ogni passo k, della direzione p(k) come l’antigradiente della F calcolato nell’iterato
k-1 -esimo, ovvero
p ( k ) = −∇F ( x (k −1) ) = − Ax ( k −1) + b = − r ( k −1) .
(5)
Questo significa che ad ogni passo il vettore p(k) coincide con la direzione di massima
pendenza per F.
In questo caso la (6) diventa
tmin
( r ( k −1) , p ( k ) )
(r ( k −1) , r ( k −1) )
=−
=
(k )
( k)
( Ap , p ) ( Ar ( k −1) , r ( k −1) )
e
x (k ) = x ( k −1) − t min r ( k −1) .
6
5. Metodo del Gradiente Coniugato
In questo metodo la scelta della direzione di discesa p(k) tiene conto non solo del
gradiente della F(x(k-1)) cioè di r(k-1) , ma anche della direzione di discesa p(k-1). In
particolare nel metodo del Gradiente Coniugato, al generico passo k, partendo dal
punto x(k-1) che è stato ottenuto muovendosi lungo la direzione p(k-1), e in cui è stato
calcolato il residuo r(k-1) (ortogonale a p(k-1)), si sceglie la nuova direzione di discesa
come quella appartenente al piano πk-1 passante per x(k-1) e individuato dai due vettori
ortogonali r(k-1) e p(k-1).
Più precisamente si ha
p ( k ) = −r ( k −1) + γ k −1 p ( k −1) , k = 1, 2,....
(6)
Il piano πk-1 interseca la funzione F(x) con un ellisse passante per x(k-1). Poiché x(k-1)
era stato scelto come punto di minimo nella direzione p(k-1) , in x(k-1) l’ellisse risulta
tangente alla direzione di discesa p(k-1) .
Poiché il punto di minimo nel piano πk-1 coincide con il centro dell’ellisse, il
parametro γk-1 sarà scelto in modo che la direzione p(k) punti verso il centro
dell’ellisse, cioè sia il coniugato, rispetto all’ellisse, di p(k-1) . Ciò significa che deve
soddisfare la seguente relazione:
( Ap( k ) , p ( k −1) ) = ( p ( k ) , Ap( k −1) ) = 0.
(7)
(k )
( k −1)
+ γ k −1 p ( k −1) nella (7) si ottiene
Sostituendo l’espressione p = − r
γ k −1
(r (k −1) , Ap (k −1) )
= (k −1)
.
(p
, Ap (k −1) )
(* *)
Utilizzando tale valore nella (6) si ottiene la nuova direzione p(k) e il nuovo punto x(k)
viene calcolato come punto di minimo nella direzione p(k), cioè
x ( k ) = x ( k −1) + λ k p (k )
(8)
con
7
(r ( k −1) , p (k ) )
λk = −
( Ap (k ) , p ( k ) )
k = 1,2,.. .
(9)
Si noti che λk in (8) gioca il ruolo di t nelle formule del paragrafo precedente e quindi
la (9) rappresenta il valore che assume λk nel punto di minimo nella direzione p(k).
Le relazioni (6), (**), (8) e (9) definiscono sostanzialmente il metodo del gradiente
coniugato.
Ci avvaliamo ora di alcuni risultati che permettono di semplificarne le formule
riducendone la complessità computazionale.
La prima semplificazione si ottiene osservando che per il residuo è possibile definire
una formula ricorsiva che lo aggiorna utilizzando una quantità che è necessaria anche
per calcolare altre grandezze, cioè
r ( k ) = Ax ( k ) − b = Ax ( k −1) − b + λk Ap ( k )
cioè
r ( k ) = r ( k −1) + λk Ap ( k ) .
(10)
La seconda semplificazione segue dal seguente
Teorema:
Nel metodo del gradiente coniugato le direzioni di discesa p(k) , con k=1,2,.., formano
un sistema di direzioni coniugate, mentre i vettori residui r(k) , con k=0,1,.., formano
un sistema ortogonale, cioè
(r(k) ,r(k-1) )=0
k=1,2,.., .
(11)
8
Osservazione:
Poiché in Rn non si possono avere più di n vettori che costituiscono un sistema
ortogonale, in linea teorica questa classe di metodi appartiene ai metodi diretti
poiché viene costruita una successione {x(k)}k=0,1,.. di vettori tali che
x ( k ) = x* = A −1b
quando k = n-1.
In pratica, però, a causa degli errori di arrotondamento, il metodo non termina al
passo k=n-1 e viene quindi utilizzato come metodo iterativo.
In molti casi comunque si verifica che il numero di iterazioni che occorrono per
raggiungere la precisione richiesta è di gran lunga inferiore alla dimensione del
sistema e questo rende il metodo molto utile per problemi di grosse dimensioni.
Utilizziamo ora il fatto che i residui in due passi successivi sono ortogonali per
ottenere un’ulteriore proprietà di ortogonalità. Infatti sostituendo nella relazione che
esprime il risultato generale,
( r ( k ) , p (k ) ) = 0
(12)
l’espressione di p(k) data dalla (11) e utilizzando la proprietà di ortogonalità (8) si
ottiene
0 = −(r ( k ) , r (k −1) ) + γ k −1 (r ( k ) , p (k −1) )
(r ( k ) , p ( k −1) ) = 0 .
Mediante questi risultati è possibile trovare una nuova espressione per il parametro λk
dato dalla (9).
Infatti poiché
(r ( k −1 ) , p ( k ) ) = − (r ( k −1 ) , r ( k −1) ) + γ k −1 ( r (k −1) , p (k −1) ) = − ( r ( k −1 ) , r (k −1) )
14
4244
3
=0
si ottiene
9
(r ( k −1) , r ( k −1) )
λk =
( Ap ( k ) , p ( k ) )
k = 1,2,..
(13)
Da questa formula si vede che se il residuo non è nullo λk è sempre positivo.
Utilizzando ora la formula ricorrente (10) per il residuo e la nuova espressione (13) di
λk è possibile trovare una formula computazionalmente più efficiente per γk-1.
Dalla formula ricorrente del residuo
r ( k −1) = r ( k −2 ) + λk −1 Ap ( k −1)
si ottiene
Ap ( k −1) =
1
( r ( k −1) − r ( k − 2) )
λk −1
e quindi
( r ( k −1) , Ap ( k −1) ) =
1
1
[(r ( k −1) , r ( k −1) ) − ( r ( k −1) , r ( k −2 ) )] =
[( r ( k −1) , r ( k −1) )] .
λk −1
λk−1
Sostituendo ora la formula di λk si ottiene
(r
( k −1)
, Ap
( k −1)
( Ap ( k −1) , p ( k −1) ) ( k −1) ( k −1)
)=
[( r
,r
)]
( k −2 )
( k −2 )
(r
,r
)
e quindi l’espressione di γk-1 diviene
? k −1
(r ( k −1) , r ( k −1) )
= ( k − 2) ( k − 2) .
(r
,r
)
(14)
10
In definitiva l’algoritmo del Gradiente Coniugato può essere schematizzato come
segue:
Scelto x(0) arbitrario, si calcola r (0)=Ax (0) - b, si prende
p(1)=- r(0)
Per k=1,2,…
(r ( k −1) , r ( k −1) )
λk =
( Ap ( k ) , p ( k ) )
x ( k ) = x ( k −1) + λk p ( k )
r ( k ) = r ( k −1) + λk Ap ( k )
If
|| r ( k ) ||22 ≥ ε
%test di convergenza non soddisfatto
% Calcolo la nuova direzione di discesa
(r ( k ) , r ( k ) )
γ k = ( k −1) ( k −1)
(r , r )
p ( k +1) = −r ( k ) + γ k p ( k )
else
%ESCI:test di convergenza soddisfatto
return x(k)
Osservazione: l’algoritmo del gradiente coniugato così ottimizzato necessita di
un’unica moltiplicazione matrice per vettore per ogni iterazione.
Velocità di Convergenza
La velocità di convergenza di un metodo iterativo si può misurare considerando di
quanto si è ridotto l’errore iniziale alla k-esima iterazione .
Per misurare l’errore si definisce la norma indotta dalla matrice simmetrica definita
positiva A su x come
11
|| x ||2A = x T Ax .
Per il metodo del gradiente Steepest Descent vale la seguente relazione
k
|| x
(k )
 K ( A) − 1 
− x || A ≤ 
 ) ≤|| x ( 0) − x * || A
 K ( A) + 1 
*
Pertanto, definendo l’errore al passo k
e (Ak ) =|| x ( k ) − x * || A
si ha
k
e(Ak )
 K ( A) − 1 
≤ 
 ) ≤ e(A0 )
 K ( A) + 1 
dove K(A) è l’indice di condizionamento di A, dato da K(A)=||A||·||A-1 ||=λmax/λmin.
 K ( A) − 1 
 ≈ 1 e quindi tanto più è lenta la
K
(
A
)
+
1


Tanto più K(A) è alto tanto più il rapporto 
convergenza.
Per il metodo del gradiente coniugato vale la seguente relazione
k
|| x ( k )
 K ( A) − 1 
 ) ≤|| x ( 0) − x * || A
− x* || A ≤ 

 K ( A) + 1 
cioè
k
 K ( A) − 1 
 ) ≤ e0A
ekA ≤ 

 K ( A) + 1 
che mostra come la convergenza di questo metodo, pur rimanendo sempre legata
all’indice di condizionamento di A sia più veloce di quella del metodo di Steepest
Descent a parità di valori di K(A).
Comunque, se la matrice A è molto mal condizionata può accadere che siano
necessari molti passi di iterazione per ottenere la convergenza.
12
Poiché l’obiettivo dei metodi iterativi è quello di ottenere una buona approssimazione
della soluzione del sistema Ax = b con, mediamente, poche iterazioni, sono state
studiate tecniche di precondizionamento che trasformano il problema originale in un
problema equivalente ma meglio condizionato.
Osservazione:
Poiché la funzione quadratica F(x) data dalla (2) assegnata la F(x)=cost rappresenta
l’espressione di un iperellissoide con eccentricità legata dal rapporto
λ max
, possiamo
λ min
dire che ad una matrice A mal condizionata corrisponde un’ iperellissoide molto
allungato, mentre ad un K(A) piccolo corrisponde un iperellissoide più arrotondato.
Metodo del Gradiente Coniugato Precondizionato
Si consideri la trasformazione
Aˆ = SAS T
detta Trasformazione di Congruenza con S non singolare e scelta in modo tale che
K ( Aˆ ) < K ( A) .
Si osserva che le trasformazioni di congruenza in genere non mantengono gli auto
valori ma ne mantengono il segno.
Allora il sistema originale Ax = b è equivalente al seguente sistema precondizionato
Aˆ xˆ = bˆ ,
con
xˆ = S −T x e bˆ = Sb .
Per questo nuovo sistema l’algoritmo del gradiente coniugato diventa
13
Scelto x̂ ( 0) arbitrario, si calcola rˆ (0) = Axˆ ( 0) − bˆ , si prende pˆ (1) = −rˆ ( 0)
Per k=1,2,…
( k −1 )
( k −1)
ˆ
ˆ
(
r
,
r
)
λˆk =
( Aˆ pˆ (k ) , pˆ ( k ) )
xˆ ( k ) = xˆ (k −1) + λˆk pˆ (k )
rˆ ( k ) = rˆ ( k −1) + λˆk Aˆ pˆ ( k )
If
|| rˆ ( k ) ||22 ≥ ε
%test di convergenza non soddisfatto
% Calcolo la nuova direzione di discesa
( rˆ ( k ) , rˆ ( k ) )
γˆ k = ( k −1) ( k −1)
( rˆ , rˆ )
pˆ ( k +1) = − rˆ ( k ) + γˆk pˆ ( k )
else
%ESCI:test di convergenza soddisfatto
return x(k)
Partendo da
x( k ) = S T xˆ ( k )
ci si chiede come è possibile esprimere il nuovo
algoritmo in termini delle grandezze originali.
Il desiderio di esprimere l’algoritmo in termini della matrice di partenza è dovuto
essenzialmente al fatto che la matrice precondizionata  presenta una struttura più
complicata e in genere risulta più densa rispetto alla matrice A. Inoltre si cerca di
risparmiare il costo della trasformazione di congruenza
Aˆ = SAS T .
14
Pertanto si sostituiscono nell’algoritmo dato le nuove grandezze con la loro
definizione:
Poiché
xˆ (k ) = S −T x (k ) , bˆ = Sb
si ha che
rˆ ( k ) = Axˆ ( k ) − bˆ = SAS T S −T x ( k ) − Sb = S ( Ax ( k ) − b ) = Sr ( k ) .
Considerando poi
xˆ ( k ) = S −T x ( k ) = S −T x ( k −1) + λk S − T p ( k ) = xˆ ( k −1) + λˆk pˆ ( k )
Si ottiene
pˆ (k ) = S −T p ( k ) .
(15)
( rˆ ( k −1) , rˆ ( k −1) )
ˆ
Per calcolare il valore λ k =
bisogna calcolare i due prodotti scalari.
( Aˆ pˆ ( k ) , pˆ ( k ) )
T
( Aˆ pˆ ( k ) , pˆ ( k ) ) = pˆ (k ) Aˆ pˆ ( k ) = ( S − T p( k ) )T SAS T S − T p ( k ) =
p( k ) S −1SAp ( k ) = p ( k ) Ap ( k ) = ( Ap ( k ) , p ( k ) )
T
T
(rˆ( k−1) , rˆ( k −1) ) = ( Sr ( k −1) )T Sr ( k −1) = r ( k−1) S T Sr ( k −1) = ( S T Sr ( k −1) )T r ( k −1)
T
da cui ponendo
~
r ( k −1) = S T Sr ( k −1)
(16)
si ha
(rˆ( k −1) , rˆ( k −1) ) = ( ~
r ( k −1) , r ( k −1) ) .
(17)
Quindi
15
( k −1)
( k −1)
( k −1)
( k −1)
~
ˆ
ˆ
(
r
,
r
)
(
r
,
r
) ~
ˆ
λk =
=
= λk .
(k )
(k )
(k )
(k )
ˆ
( Apˆ , pˆ ) ( Ap , p )
Ora è possibile scrivere il nuovo iterato del metodo precondizionato come
~
xˆ ( k +1) = xˆ ( k ) + λk pˆ ( k )
da cui moltiplicando per ST , si ha
~
x ( k +1) = x ( k ) + λk p ( k ) .
Analogamente considerando il residuo si ha
~
rˆ( k +1) = rˆ (k ) + λk +1 Aˆ pˆ ( k +1)
moltiplicando per S -1 e ricordando che rˆ
(k )
= S r (k ) si ottiene
~
r ( k +1) = r (k ) + λk +1S −1SAS T S −T p (k +1)
cioè
~
r ( k +1) = r (k ) + λk +1 Ap( k +1) .
Per calcolare γ k utilizziamo la (17) e otteniamo
(~
r (k ) , r (k ) )
~
γ k = ~ ( k −1) ( k −1) .
(r , r )
Infine consideriamo
pˆ ( k +1) = −rˆ ( k ) + γ~k pˆ ( k )
e moltiplicando per S T si ha
S T S −T p ( k +1) = − S T Sr ( k ) + γ~k S T S −T p ( k )
da cui segue che
p ( k +1) = −~
r ( k ) + γ~k p ( k ) .
16
Tutti questi calcoli ci hanno permesso di ottenere le formule di base per il metodo del
gradiente coniugato precondizionato, scritto in termini delle grandezze originali più
un nuovo vettore
~
r = S T Sr .
Se conoscessimo esplicitamente le matrici S e ST che compaiono nella trasformazione
di congruenza
Aˆ = SAS T sarebbe facile calcolare ~
r.
In realtà le matrici S e ST sono definite implicitamente, in quanto assegnata A
simmetrica e definita positiva, per ottenere un buon precondizionatore, noi cerchiamo
una matrice M simmetrica definita positiva che approssimi A, tale che
M-1A ≅ I.
Poiché M è simmetrica e definita positiva essa può essere definita come
( S T S ) −1 = M
in quanto S è non singolare.
Perciò si ottiene che
~
r ( k ) = S T Sr ( k ) = M −1r ( k )
cioè
M~
r (k ) = r ( k ) .
Questo significa che per calcolare il nuovo residuo occorre risolvere questo sistema
lineare con la matrice M che approssima A.
L’algoritmo del gradiente coniugato precondizionato diviene quindi:
17
Scelto
x ( 0)
arbitrario,
si
calcola
r (0) = Axˆ ( 0) − bˆ ,
si
calcola
M~
r ( 0) = r ( 0) e si prende p (1) = −~
r ( 0)
Per k=1,2,…
~ (~
r ( k −1) , r ( k −1) )
λk =
( Ap ( k ) , p ( k ) )
~
x ( k ) = x ( k −1) + λk p ( k )
~
r ( k ) = r ( k −1 ) + λk Ap ( k )
Risolvi
If
M~
r ( k ) = r (k )
(~
r (k ) , r (k ) ) ≥ ε
%test di convergenza non soddisfatto
% Calcolo la nuova direzione di discesa
(~
r (k ) , r (k ) )
~
γ k −1 = ~ ( k −1) ( k −1)
(r
,r
)
p ( k +1) = − ~
r ( k ) + γ~ p ( k )
k
else
%ESCI:test di convergenza soddisfatto
return x(k)
Si osserva che la differenza sostanziale dell’algoritmo precondizionato sta nella
risoluzione, ad ogni passo, di un sistema lineare; è evidente che per M=I l’algoritmo
si riduce a quello non precondizionato.
Ora quindi il problema diviene quello di costruire una matrice M che approssima A e
che sia facile da invertire.
Prima però vogliamo ci soffermiamo a capire perché dal problema precondizionato
Aˆ xˆ = bˆ
con
Aˆ = SAS T
18
siamo passati al problema di costruire la matrice
M = ( S T S ) −1
con M matrice che approssima A.
Â
Questo è giustificato dal fatto che, se considero di premoltiplicare
per
S −T
(avendo supposto S invertibile) ottengo
S T Aˆ S −T = S T SA = M −1 A
ˆ è simile a M −1 A e quindi la trasformazione di congruenza S sarà tanto
cioè A
migliore quanto più riduce l’indice di condizionamento di A, cioè
λ ( M −1 A)
cond ( Aˆ ) = cond ( M −1 A) = max
.
λmin ( M −1 A)
−1
Quindi tanto più M approssima A, tanto più c’è la speranza che M A approssimi
l’identità, cioè sia con un indice di condizionamento prossimo ad 1.
19
Calcolo di possibili precondizionatori
In letteratura esistono numerose possibilità di precondizionatori e il problema è
ancora oggetto di ampie ricerche.
1. Precondizionatori Diagonali
Si definisce la matrice M come la matrice diagonale con gli elementi di A sulla
diagonale, cioè
M = Diag (a11, ......., a nn ) = D A
e quindi
M −1 = Diag (
1
1
.,......,
) = D A−1 .
a11
a nn
Questa scelta corrisponde a considerare la matrice S che compare in
Aˆ = SAS T
T
−1
e che corrisponde a M = ( S S ) pari a
S = S T = Diag (
1
.....
a11
1
−1
) = DA 2 .
a nn
r = r facilmente
Questa scelta è molto semplice in quanto rende il sistema M ~
risolvibile; purtroppo però non sempre porta a buoni risultati.
2. Precondizionatori basati sulla fattorizzazione di Cholesky incompleta
Poiché la matrice A è simmetrica e definita positiva, l’idea è di considerare una
matrice triangolare inferiore G che presenti una qualche struttura sparsa e
approssimi al meglio possibile il fattore triangolare inferiore L della
fattorizzazione di Cholesky di A=LLT e scegliere
M=GGT.
20
Questa scelta ci assicura che M è simmetrica definita positiva e il sistema
M~
r = r è facilmente risolvibile data la fattorizzazione M=GGT e la sparsità dei
fattori triangolari G e GT.
La matrice S che corrisponde a questa scelta di M risulta essere
S = G −1 .
Infatti
M −1 = G −T G −1 = S T S
quanto più G ≅ I.
Il calcolo di G viene fatto in genere eseguendo una fattorizzazione parziale di A
che consiste nel porre a 0 i valori di G che corrispondono ai valori aij nulli e
calcolando solo gli altri.
Osservazione finale:
Il metodo del gradiente richiede che la matrice A sia simmetrica e definita positiva.
Esiste però la possibilità di estenderlo a matrici qualsiasi purchè non singolari.
Dato il sistema lineare
Ax = b
lo si sostituisce con il sistema
ATAx = ATb
In cui la matrice ATA è simmetrica e definita positiva. Ora è possibile applicare il
metodo del gradiente.
Per valutare l’algoritmo occorre però considerare sia il costo superiore in termini di
complessità computazionale, dovuto all’introduzione dei prodotti matrice per vettore,
sia la convergenza più lenta poichè l’indice di condizionamento è più elevato dovuto
al fatto che K2(ATA) = (K2(A))2.
21