Matematica Computazionale(6cfu) Ottimizzazione · 2016. 5. 13. · Matematica Computazionale,...

1

Corso di Laurea in Infomatica

Corso di Laurea in Matematica

Matematica Computazionale(6cfu)

Ottimizzazione(8cfu)

(a.a. 2015-16, lez.13)

Docente: Marco Gaviano

(e-mail:[email protected])

2

Matematica Computazionale, Ottimizzazione, a.a. 2015-16, Lezione, n.13

Minimizzazione unidimensionale

Un problema di minimizzazione in una variabile si può formulare come

minimizzare f() R (1)

E' importante poter disporre di buoni algoritmi che risolvono (1)

efficientemente . Ciò perché la risoluzione in due o più dimensioni si

baserà su minimizzazioni unidimensionali.

Ipotesi 1: f() abbia derivata prime e seconde nell'intervallo [a,b]; sia

inoltre strettamente convessa, ovvero f =f"()>0.

Esempio di funzione

strettamente convessa,

3

Matematica Computazionale, Ottimizzazione , a.a. 2015-16, Lezione, n.13

Si danno ora alcuni algoritmi per la sua risoluzione.

Algoritmo della sezione aurea (golden section)

Si parte da un intervallo [ai.bi] in cui esiste un punto di minimo ed in esso

si scelgono i punti v e w.

Se f(v)f(w), allora il minimo si trova nell’intervallo [ai+1,bi+1] [ai w]; se

f(v)>f(w), allora il minimo si trova nell’intervallo [ai+1,bi+1] [v bi] (La

figura illustra il procedimento).

Si ripete la procedura a partire da [ai+1 bi+1] e si ottiene [ai+2 bi+2].

Si ripete l’iterazione finché la misura dell’intervallo è sufficientemente

piccola.

4

Matematica Computazionale, Ottimizzazione , a.a. 2015-16, Lezione, n.13

bi v w ai

v w ai+1 bi+1

caso f(v) > f(w)

w bi+1 v

v w ai bi

ai+1

caso f(v)f(w),

5


Algoritmo 1 (sezione aurea) >0. F1 =(3-51/2)/2, F2 =(51/2 -1)/2,

Sia [a,b]= [a0 ,b0 ], i=0,

v= ai+F1(bi -ai ), w= ai+F2(bi -ai )

fv=f(v), fw=f(w),

while |(bi -ai )|>

if fv fw

ai+1 = ai ,bi+1 =w,

w=v, v= ai+1+ F1 (bi+1 –ai+1), fw=fv, fv=f(v) ,

else

ai+1 = vi ,bi+1 =bi,

v=w, w= ai+1+ F2 (bi+1 –ai+1), fv=fw, fw=f(w)

end

i=i+1

end

= (ai + bi)/2, (approssimaziome punto di minimo)

6


Convergenza

Nel caso si ponga =0 allora l’algoritmo della sezione aurea genera una

successione infinita {[ai ,bi ]}. Può allora dimostrarsi la seguente,

Proposizione

L'ipotesi 1 sia valida. Allora

con =(51/2 -1)/2=0.6180,

ia)θ(b2

15a)(bab

i

ii

7


Algoritmo dell’interpolazione parabolica

Si costruisce una successione i . Si parte

dal valore della f() nei punti i e i-1 e dal

valore della sua derivata nel punto i-1. Si

costruisce la parabola (polinomio di 2° grado)

che nei punti i e i-1 ha lo stesso valore di

f() ed inoltre in i-1 ha la stessa derivata. Il

minimo della parabola è utilizzato per

costruire il nuovo punto i+1.

i i-1 i+1

parabola

f()

minimo parabola

8


Algoritmo 2 (interpolazione parabolica) >0. f '(0 =0)<0, 1 >0,

Poni i=1,

while |f (i)|> ,

end

=i+1 (approssimaziome punto di minimo)

))f(λ)(f(λ)(λ' )fλ(λ

)λ)(λ(λ' fλ

2

1λ

i1ii1ii

2

1ii1i1i1i

9


Convergenza

Come nel caso dell’algoritmo della sezione aurea, l’algoritmo

dell’interpolazione parabolica genera una successione infinita {i}.

Può allora dimostrarsi la seguente

Proposizione

L'ipotesi 1 sia valida. Allora la successione {i } converge al punto di

minimo * di f(), inoltre ,

con = 20.5 e (0,1) ed opportune costanti.

iγ

i ωθ |λ*λ|

10


Algoritmo dell’interpolazione cubica

Si costruisce una successione i . Si parte

dal valore della f() e della sua derivata nei

punti i e i-1 . Si costruisce la cubica

(polinomio di 3° grado) che nei punti i e i-1

ha lo stesso valore di f () ed della sua

derivata. Il minimo della cubica è utilizzato

per costruire il nuovo punto i+1.

i i-1 i+1

cubica

f()

minimo cubica

11


Algoritmo 3 (interpolazione cubica) >0. f '(0 =0)<0, 1 >0,

Poni i=1,

while |f '(i)|>

end

=i+1 (approssimaziome punto di minimo)

0.5

1ii

2

11ii2

i1-i1iii1i1

21ii

21i1iii1i

))(λ' )f(λ' f)(uλsgn(λu

)λ))/(λλ3(λ-)(λ' f)(λ' fu

con

u)(λ' f)(λ' f

)uu)(λ')(fλ(λλλ

12


Convergenza

L’algoritmo dell’interpolazione parabolica genera una successione

infinita {i}. Può allora dimostrarsi la seguente

Proposizione

L'ipotesi 1 sia valida. Allora la successione {i } converge al punto di

minimo * di f(), inoltre

con = 2 e (0,1) ed opportune costanti.

iγ

i ωθ|λ*λ|

13


Minimizzazione in n dimensioni,

Come avviene spesso nell'Analisi Numerica si cerca di raggruppare gli

algoritmi che risolvono un certo problema in classi o famiglie.Vengono

quindi studiate le proprietà (convergenza, tasso di convergenza, ...) della

classe. La maggior parte degli algoritmi di minimizzazione per il

problema

min f(x) xRn (2)

possono inquadrarsi nel seguente schema generale o modello di

algoritmo. Da notare che gli algoritmi sono iterativi e la soluzione

viene trovata come punto limite di una successione {xk}

14


Algoritmo 4 (schema,n-dimensional minimization). >0

Scegli x0 Rn, x0 starting point; k=0;

while ||f(xk)|| >

calcola pkRn , pk la direzione di ricerca tale che

-(pk)Tf(xk ) > 0(direzione di discesa);

calcola un numero positivo k, (la lunghezza del passo),

(step length), tale che,

f(xk +kpk ) < f(xk );

calcola il nuovo punto

xk+1=xk+kpk

k=k+1

end

x*=xk, (approssimaziome punto di minimo).

15

L'algoritmo 4 è abbastanza generale e comprende la maggior parte degli

algoritmi di minimizzazione. Il suo comportamento è illustrato nella

seguente figura.


-f(xk) xk pk

pkdirezione di ricerca

f(xk)=gradiente di f(x)

partendo da xk si procede

lungo pk scegliendo un

passo k, e xk+1= xk+ k pk

e f(xk+1) f(xk). xk+1

16


Sostanzialmente da un punto xk si procede lungo una direzione di

discesa per la f(x). Lungo pk si sceglie un passo opportuno. Si calcola

xk+1 e si ripete la procedura fino a che il criterio di stop è soddisfatto.

Ciò che distingue i vari algoritmi è la scelta e della direzione di ricerca

pk e quella del passo k. Il criterio di stop usato nella struttura while è,

|f(xk )| (3)

con un numero positivo prefissato; ma può usarsi un altro criterio di

stop.

17


La scelta di k

Il valore ottimale k è quello per cui vale la seguente relazione,

f(xk +kpk) = min f(xk +pk), 0. (4)

Si parla allora di ricerca lineare esatta. Chiaramente questa scelta

assicura il massimo decremento della funzione lungo pk. Per il calcolo di

k si potrebbe dunque pensare di utilizzare l’algoritmo della sezione

aurea o dell’interpolazione parabolica o cubica. In tal caso la funzione

unidimensionale cambia ad ogni iterazione, precisamente si ha

f k() = f(xk +pk),

La scelta di un 'piccolo' negli algoritmi unidimensionali garantisce

che (4) è soddisfatta con una buona approssimazione.

18


Questa scelta, benché corretta, è inefficiente nel senso che essa richiede

un numero molto alto di valutazioni di funzione (Un criterio di

valutazione della bontà di un algoritmo è il numero di valutazioni di

funzione che la sua applicazione richiede; nel Calcolo Numerico questo

criterio è usato spesso). Insistere sulla ricerca lineare in modo da

ottenere grande accuratezza porta ad un risultato disastroso.

Gli algoritmi 1,2 e 3, pertanto, vengono utilizzati con criteri di stop

basati sul principio di Goldstein-Armijo. Vale a dire si sceglie un k tra

tutti gli che soddisfano la seguente relazione,

0< -1(pk)Tf(xk) f(xk) - f(xk +pk) -2(pk)

Tf(xk) (5),

con 1 e 2 numeri positivi tali che 0<1 2<1,

19


La figura illustra il significato di (5).

2(pk)Tf(xk) f(xk +pk)- f(xk) 1(pk)

Tf(xk)

1(pk)Tf(xk)

f(xk +pk)- f(xk)) 2(pk)

Tf(xk)

Il principio di Goldstein-Armijo

20


Il caso 1 = 2= 0.5 equivale come costo computazionale ad effettuare

la ricerca lineare esatta ossia a calcolare k secondo (4).

Un altro criterio è il seguente, si sceglie k in modo da soddisfare, la

relazione,

|(pk)Tf(xk +pk)|(pk)

Tf(xk) (6)

con 01. Questa è giustificata dal fatto che nel punto di minimo

lungo pk

(pk)T f(xk +pk)=0

Dunque si rende meno stringente la condizione di convergenza sul passo.

21


La scelta della direzione di ricerca pk

Per la scelta della direzione di ricerca pk nell’algoritmo 4 sono stati

proposti centinaia di algoritmi che possono raggrupparsi in classi. Si

considererà il metodo del gradiente, il metodo di Newton, i metodi quasi-

Newton ed il metodo del gradiente coniugato. La proprietà che si

richiede a pk è che sia di discesa; vale a dire,- (pk)T f(xk)>0. Questa

proprietà è abbastanza debole e da sola, anche con una scelta ottimale del

passo, non può garantire la convergenza dell’algoritmo 4. La condizione

che quindi spesso si richiede è

-(pk)T f(xk )>||pk|| || f(xk )|| (7)

con (0,1). Questa assicura che l'angolo tra il gradiente negativo e pk

non tenda a /2 (vedi figura).

22


-f(xk) xk

pk

pkdirezione di ricerca

f(xk)=gradiente di f(x)

questo angolo non

deve tendere a /2

23


Il metodo del gradiente o di più ripida discesa

La (pk)Tf(xk ) con ||pk||=1 dà la derivata prima di f(xk+pk ) (vista

come funzione di ) per =0. E' dunque naturale scegliere come

direzione di ricerca quella direzione che rende minima pTf(xk)/||p|| .

Si consideri il problema

in cui la norma è quella Euclidea, (||x||=(x1 ,x2 ,...,xn)0.5 ), allora la

soluzione di g(p)=0 è

p=-f(xk )

La scelta pk=-f(xk ) nell’algoritmo 4 dà quindi origine al metodo di più

ripida discesa.

)f(xp

pg(p) k

T

Rpmin

n

24


Vale la seguente proposizione

Proposizione

f(x) sia una funzione con derivate prime e seconde continue. Il metodo

del gradiente con una opportuna scelta del passo (utilizzo teoricamente

del passo ottimale o più realisticamente degli algoritmi 1, 2, e 3) generi

una successione convergente a x* . Allora f(x*)=0 ovvero x* è un

punto di stazionarietà.

Se vale la condizione di uniforme convessità

m||x||2xTH (x)x M||x||2, xRn

con m ed M, mM, numeri reali positivi e H(x), l’hessiano di f(x) calcolato in x, si ha convergenza lineare,

||xk-x*||ck

per k=1,2,... , c un opportuno numero reale positivo e (1-m/M)0.5.

25


Il metodo di più ripida discesa per una funzione di due variabili,

nel caso il punto di stazionarietà sia minimo locale

-f(•)

xk

xk+1

26


Comportamento del metodo del gradiente per m<<M e mM.

Nel secondo caso il problema può dirsi malcondizionato e risulta

evidente l'inefficienza del metodo del gradiente.

mM

m<<M

27


Il metodo di Newton

Il metodo di Newton per problemi di minimo è derivato dallo stesso

metodo per la risoluzione di sistemi di equazioni nonlineari. Esso nel

caso sia applicabile e convergente ha un tasso di convergenza quadratico.

Pertanto in poche iterazioni permette di avere una buona

approssimazione della soluzione. Ma essendo computazionalmente

costoso esso non viene applicato spesso, è comunque diventato il

modello per altri metodi che hanno un tasso di convergenza comparabile

ma che richiedono ad ogni iterazione una mole di calcolo notevolmente

inferiore. La giustificazione del metodo è la seguente.

Si consideri il seguente sviluppo di Taylor di f(x)

f(xk +p)= f(xk )+p Tf(xk ) +½pTHkp+...

con Hk=H(xk), H(x) l'Hessiano di f(x) calcolato in x .

28


Si approssimi ora la f(x) con la somma dei primi tre termini e si

minimizzi rispetto a p. Si trova che p deve soddisfare la relazione

Hkp = -f(xk ).

Questa suggerisce di calcolare pk nell’algoritmo 4 secondo la formula,

pk = -(Hk)-1f(xk ) (8)

Il metodo che ne deriva è chiamato metodo di Newton. La scelta di p

tiene conto di informazioni sulla derivata seconda di f(x).

La scelta pk = -f(xk ) nel metodo del gradiente era basata solo su

informazioni del primo ordine. Dunque con la scelta (8) è naturale avere

risultati migliori per ciò che riguarda il tasso di convergenza.

29


Si possono fare le seguenti osservazioni.

• Nel caso di funzioni quadratiche con k =1 una sola iterazione è

sufficiente ad individuare il minimo.

• Il costo numerico per iterazione è molto alto in quanto si richiede

e il calcolo dell'hessiano e il calcolo della sua inversa.

• La direzione p k è di discesa solo se la matrice H k è definita

positiva per ogni k.

• La condizione (7) risulta soddisfatta solo se H k per k, non

tende a diventare semidefinita positiva ovvero vale la condizione

di uniforme convessità.

m||p||2pTH k p M||p||2, pRn

per k=1,2,... e con m ed M, mM, numeri reali positivi.

30

Per l'algoritmo 4 con la direzione di ricerca data da (8) si ha la seguente

proposizione.

Proposizione

f(x) sia una funzione con derivate prime e seconde continue ed inoltre

m||x||2xTH(x) x M||x||2 xRn

con m ed M, mM,numeri reali positivi, allora il metodo di Newton con

una scelta opportuna del passo (utilizzo del passo ottimale o più

realisticamente degli algoritmi 1, 2 e 3) genera una successione

convergente al minimo x* . Si ha inoltre convergenza quadratica,

k=1,2,..., e c un opportuno numero reale positivo e (0,1).


k

cxxk

2*

31


Si noti che nel metodo di Newton il passo ottimale (4) k tende a 1 al tendere di xk alla soluzione.

La proposizione mette in evidenza le restrizioni che devono valere per

l'applicazione del metodo di Newton. In realtà le ipotesi fatte possono

restringersi ad un intorno di un punto di minimo locale, la convergenza è

comunque dimostrabile solo partendo da un punto in tale intorno.

Nel caso di funzioni generali in cui o l'hessiano è definito negativo o la

matrice H(x) è singolare, allora chiaramente (8) o non è applicabile o

genera un direzione non di discesa. In tal caso sono stati proposti nella

letteratura numerosi metodi chiamati di Newton modificato, nei quali la

pk è sempre di discesa. Generalmente se si è lontani dalla soluzione si

applica Newton modificato (eventualmente si può utilizzare la direzione

del gradiente negativo), nelle vicinanze della soluzione si passa alla (8).

32


Un esempio di calcolo della direzione di ricerca è dato da

if f(xk )T (Hk)

-1f(xk ) ||f(xk )|| ||(Hk)-1f(xk )|| oppure

l’hessiano è non-singolare e definito positivo allora,

pk =-(Hk)-1f(xk )

else

pk =- f(xk )

end.

con (0,1).

33


Metodi quasi Newton

I metodi quasi Newton, o a metrica variabile come anche vengono

chiamati, sono stati proposti per ovviare agli inconvenienti del metodo di

Newton già citati, ovvero il calcolo dell'hessiano e della sua inversa. In

questi metodi non si vuole utilizzare l'hessiano ma solo i valori della

funzione f(x) e i valori del gradiente f(xk) già calcolati. Man mano che

si procede con le iterazioni si cercherà di approssimare mediante una

matrice Bk l’inverso dell'hessiano H(xk). Vale a dire si cerca di

soddisfare il limite

limk Bk=H(x*)-1.

dove con x* si indica il minimo a cui la successione {xk} tende. La

direzione di ricerca si calcola quindi secondo una formula uguale alla (8)

ovvero

pk =-Bkf(xk )

34


Come può costruirsi la matrice B?

Nel caso f(x) sia una funzione quadratica

f(x) = dTx+(xTHx)/2 (9)

con H matrice nn simmetrica definita positiva e d un vettore in Rn allora vale la seguente relazione

Hsk=f(xk+sk)- f(xk)

ovvero

sk= H-1(f(xk+sk)- f(xk)). (10)

Per una funzione non quadratica la (10) non vale. Ma se si è in un

intorno sufficiente piccolo di un punto di minimo locale con hessiano

definito positivo si ha

sk H(xk )-1(f(xk+sk)- f(xk)).

dove si può porre sk=xk+i-xk.

35


La (10) è dunque utilizzata per la costruzione della successione delle

matrici Bk. Si costruisce la Bk+1 a partire dalla Bk, mediante la formula

Bk+1 =Bk+Uk.

La correzione Uk è aggiunta a Bk in modo che valga la (10) ossia,

(Bk+Uk) –1(xk+1-xk)=(Bk+1)

-1 (xk+1-xk)= f(xk+1)- f(xk).

Esistono numerose espressioni per Uk. Un esempio è la seguente,

con sk = (xk+1 -xk) e yk = f(xk+1)- f(xk) e v un qualsiasi vettore tale che

vTyk0.

k

T

T

kkkk1k

yv

)vyB(sBB

36


Se si impone la simmetria su Uk allora si ha la formula di Broyden di

rango 1

Se si pone

con v e z vettori arbitrari ma tali che vTyk0, zTyk0, si ha un'altra

famiglia di formule per le quali (10) è valida. La scelta in (11) di v = sk e

z = Bk yk dà

chiamata formula di Davidon-Fletcher-Powell che è una formula

simmetrica. Inizialmente si pone B0 =I, con I la matrice identità .

(Broyden). y)yB(s

)yB)(syB(sBB

k

T

kkk

T

kkkkkkk1k

)(11 yz

zyB

yv

vsBB

k

T

T

kk

k

T

T

kkk1k

Powell)Fletcher(Davidon yBy

)y(ByB

ys

ssBB

kk

T

k

T

kkkk

k

T

T

kkk1k

37


Per questa formula si hanno i seguenti risultati.

Proposizione

L'algoritmo 4 applicato alla funzione quadratica (9) con =0 e pk

calcolato mediante

pk =-Bkf(xk ),

Bk mediante la formula di Davidon-Fletcher-Powell ed k ottimale,

converge al punto di minimo in un numero finito di passi minore od

uguale ad n.

Per funzioni non quadratiche si ha

38


Proposizione

Sia f(x) una funzione avente derivate prime e seconde continue; inoltre

per xC(x0) (di misura finita)

C(x0)= {x | f(x)f(x0) }

e per ogni yRn sia

m||y||2yTH(x)y M||y||2

con m, M costanti positive, allora l'algoritmo 4 con pk =-Bkf(xk ), Bk

dato dalla formula di Davidon-Fletcher-Powell, k ottimale e =0

partendo da x0 o si ferma dopo un numero finito di iterazioni trovando il

punto di minimo o genera una successione convergente al punto x* che

minimizza f(x) su C(x0).

39


Convergenza del metodo quasi-Newton con la formula Davidon-

Fletcher-Powell

punto di partenza x0

punto limite x*

(minimo locale)

40


Per il tasso di convergenza si ha

Proposizione

Siano valide tutte le ipotesi della proposizione precedente allora la

successione {xk} converge superlinearmente al punto di minimo x*,

ovvero,

0*

*lim

1

xx

xx

k

k

k

41


Metodo del gradiente coniugato

I metodi del gradiente coniugato furono introdotti da Fletcher e Reeves

con lo scopo di ridurre il numero di operazioni aritmetiche richiesti dai

metodi quasi-Newton. Nei primi metodi introdotti nella letteratura la

direzione di ricerca è calcolata secondo la formula

pk=-f(xk)+kpk-1 (12)

Vale a dire pk è ottenuto come combinazione lineare del gradiente

negativo e della direzione di ricerca precedente. Il parametro è scelto in

diversi modi. Si hanno le formula di Fletcher-Reeves e Polak-Ribiere

per

(13) Reeves)-(Fletcher )f(x)f(x

)f(x)f(xγ

1k

T

1k

k

T

kk

(14) Ribiere)-(Polak )f(x)f(x

)f(x))f(x)f(x(γ

1k

T

1k

k

T

1kkk

42


Vediamo ora la motivazione del metodo. Ciò giustificherà anche il nome

'gradiente coniugato'. Consideriamo la definizione di direzioni coniugate.

Definizione

Sia data una matrice nn H simmetrica e definita positiva, m vettori pi ,

i=1,2,...,m, si dicono coniugati rispetto ad H se vale la relazione,

La scelta di k in (13) e (14) è stata fatta in modo che nel caso si debba

minimizzare la funzione quadratica (9) con la scelta ottimale del passo k

allora le prime n direzioni di ricerca risultano coniugate rispetto ad H.

Può dimostrarsi

Proposizione

L'algoritmo 4 con pk calcolato secondo (12) e k dato da (13) o (14),

con scelta ottimale del passo e =0, qualora applicato ad una funzione

quadratica converge al minimo in un numero di iterazioni mn.

.1,2,...mi j,iper 0,Hpp j

T

i

43


Per l'applicazione ad una funzione non quadratica si ha

Proposizione

Sia f(x) una funzione avente derivate prime e seconde continue; inoltre

per xC(x0), con

C(x0)= {x | f(x)f(x0)}

e per ogni yRn sia,

m||y||2yTH(x)y M||y||2

con m, M costanti positive, allora l'algoritmo 4 con pk calcolato secondo

(12) e k dato da (13) o (14) e =0, partendo da xo o si ferma dopo un

numero finito di iterazioni trovando il punto di minimo o genera una

successione convergente al punto x* che minimizza f(x) su C(xo)

44


I metodi del gradiente coniugato sono usualmente applicati

riinizializzando dopo n o 2n,…, iterazioni la direzione di ricerca.

Usualmente si pone all'inizio di ogni ciclo pk=-(xk).

La riinizializzazione è motivata dal fatto che si può dimostrare che n

iterazioni di un metodo del gradiente coniugato equivalgono ad una

iterazione del metodo di Newton. I metodi del gradiente coniugato hanno

tasso di convergenza inferiore a quello dei metodi Quasi-Newton. Ma

sono semplici da applicare e non richiedono la memorizzazione di una

matrice (Quelli proposti più recentemente richiedono la memorizzazione

di più vettori), pertanto in problemi di grandi dimensione sono i metodi

più efficienti.

L'algoritmo che ne consegue è il seguente

45


Algoritmo (Polak-Ribiere , gradiente coniugato con restart). >0.

Scegli x0Rn, (starting point); k=0; p-1=(0,…,0);

while |f(xk)|>

pk=-f(xk)+kpk-1

calcola k tale che f(xk +kpk) = min f(xk +pk), 0

xk+1=xk+kpk

k=k+1,

end

x*=xk, (approssimaziome punto di minimo).

altrimenti)f(x)f(x

)f(x))f(x)f(x(

n di multiploun èk se 0

γ

1k

T

1k

k

T

1kkk

Matematica Computazionale(6cfu) Ottimizzazione · 2016. 5. 13. · Matematica Computazionale,...

Documents

Transcript of Matematica Computazionale(6cfu) Ottimizzazione · 2016. 5. 13. · Matematica Computazionale,...