+ All Categories
Home > Documents > Gradient Conjugat

Gradient Conjugat

Date post: 06-Nov-2015
Category:
Upload: iulian-boboc
View: 266 times
Download: 4 times
Share this document with a friend
Description:
Referat Analiza Numerica
27
CALCUL PARALEL METODA DE GRADIENT CONJUGAT -Tema de casa-
Transcript

Calcul Paralel - Metoda de Gradient Conjugat

CALCUL PARALELMETODA DE GRADIENT CONJUGAT

-Tema de casa-

Boboc Iulian & Andone Alexandru 2131CIntroducere

Metodele de optimizare sunt n general metode de descretere, ce determin minimul unei funcii U de n variabile reale care se numete funcie scop sau funcie obiectiv. De aici i denumirea lor de metode de minimizare a funciilor de mai multe variabile. Evident, problema gsirii maximului revine la minimizarea funciei cu semn schimbat. Metodele de descretere au convergena global, adic permit gsirea soluiei chiar dac punctul de plecare este ndeprtat de soluie.

Metodele de optimizare au un domeniu de aplicabilitate foarte larg. Pe de o parte, majoritatea fenomenelor naturii sau economice reprezint compromisuri ntre cauze contradictorii, i ca atare multe din problemele ingineriei, economiei, matematicii, statisticii, medicinei, dar mai cu seam procesele decizionale se pot formula ca probleme de optimizare. Pe de alt parte, majoritatea metodelor numerice pot fi reformulate ca probleme de optimizare. Aceste reformulri duc uneori la obinerea unor metode performante, cum ar fi cele pentru rezolvarea sistemelor de ecuaii liniare, i cele pentru rezolvarea sistemelor de ecuaii neliniare

Metodele de gradient conjugat reprezint o contribuie important n panoplia metodelor de optimizare fr restricii de mari dimensiuni. Algoritmii asociai acestor metode sunt caracterizai de cerine modeste de memorie i au proprieti foarte bune de convergen global. Popularitatea lor este datorat n parte simplitii expresiei lor algebrice, implementrii foarte rapide i uoare n programe de calcul, precum i eficienei lor n rezolvarea problemelor cu un numr mare de variabile.

Metodele de gradient conjugat au fost proiectate de Magnus Hestenes (1906-1991) i Eduard Stiefel (1909-1978) n lucrarea Methods of conjugate gradients for solving linear systems, J. Research Nat. Bur. Standards Sec. B. 48, 409-436 (1952) n care au prezentat un algoritm pentru rezolvarea sistemelor algebrice liniare, cu matricea simetric i pozitiv definit. n 1964, metoda a fost extins de Fletcher i Reeves la probleme de optimizare neliniar fr restricii. De atunci, un numr foarte mare de algoritmi de gradient conjugat au fost elaborai. Chiar dac algoritmii de gradient conjugat au peste 50 de ani de existen, totui acetia continu s fie de un considerabil interes n particular datorit proprietilor lor de convergen i eficienei n rezolvarea problemelor de optimizare de mari dimensiuni.

Metodele de gradient conjugat nu se deosebesc esenial de metodele cvasi-Newton din punct de vedere al scopului, i anume obinerea minimului unei forme ptratice n n iteraii. Ambele clase de metode necesit calculul derivatelor pariale de ordinul nti i au aceeai convergen superliniar. Deosebirea esenial const n faptul c metodele de gradient conjugat nu necesit memorarea unei matrice.

Consideraii teoretice

Metodele de gradient. Aceste metode sunt tipic de ordinul I i sunt caracterizate prin alegerea n fiecare punct curent vk a unei direcii de deplasare dk opus gradientului local:

dk = -gk

(1)

Dezvoltnd n serie Taylor funcia n vecintatea punctului vk i reinnd doar termenii de ordinul nti, se obine:

(2)

ns,

(3)

egalitatea avnd loc numai n cazul (1). Prin urmare, pentru orice pas sk >0 alegerea direciei de cutare conform relaiei (1) asigur local descreterea maxim posibil a funciei obiectiv f.

Algoritmul de calcul al punctului de minim prin metoda gradientului este prezentat n cele ce urmeaz:

0. Se alege un punct iniial astfel nct mulimea

s fie mrginit.

1. Se iniializeaz k = 0.

2. Se calculeaz gk = f v (vk ) .

Dac || g k ||= 0 , atunci i algoritmul este oprit.Altfel, se alege direcia dk = -gk i se trece la pasul 3.

3. Se determin pasul sk .

4. Se calculeaz vk+1 = vk + sk dk , se actualizeaz k prin nlocuirea lui k cu k+1 i se revine la pasul 2.

Metodele de gradient conjugat.Aceste metode au principalul avantaj c au o convergen bun, iar numrul de operaii aritmetice necesare pe iteraie este relativ redus. Prin urmare, aceste metode sunt utilizabile i n rezolvarea unor probleme de optimizri de dimensiuni mari.

Direcia de cutare pentru metodele de gradient conjugat este obinut pe baza relaiei (4):

(4)

cu d0 = - g0, iar parametrul scalar k este specific metodei i contribuie la accelerarea vitezei de convergen.

Parametrul k are urmtoarele expresii posibile care conduc la diverse variante de metode de gradient conjugat:

- pentru metoda Fletcher-Reeves:

(5)

- pentru metoda Polak-Ribire:

(6)

- pentru metoda Hestenes-Stiefel:

(7)

Din motive de convergen, pentru implementrile practice este recomandat reiniializarea algoritmului dup un numr de l n + 1 iteraii, folosindu-se l = 0. Din aceleai motive este recomandat utilizarea unei proceduri de cutare liniar exact pentru determinarea lungimii sk a pasului.

Comparaie ntre gradientul descendent (cu verde) i gradientul conjugat (n rou) pentru minimizarea funciei asociat unui sistem liniar. Gradientul conjugat convergete n maxim n pai, n exemplul de mai sus n 2 pai. (unde n este mrimea matricii de sistem)Algoritmul de calcul al punctului de minim

EMBED Equation.3 prin metoda gradientului conjugat este prezentat n cele ce urmeaz i are urmtorii pai:

0. Se alege un punct iniial

1. Se iniializeaz k = 0.

2. Se calculeaz g0 = f v (v0 ) .Dac || g 0 || = 0 , atunci i algoritmul este oprit.

Altfel, se alege direcia d0 = - g0 i se trece la pasul 3

3. Se determin pasul optimal

4. Se calculeaz v k+1 = vk + skdk

5. Se calculeaz g k+1 = f v (v k+1 ) Dac || g k+1 || = 0 , atunci i algoritmul este oprit.

Altfel, se trece la pasul 6.

6. Se calculeaz k+1 cu una din formulele cunoscute i direcia de cutare cu formula

d k+1= - g k+1 + k+1 d k

Apoi, se actualizeaz k prin nlocuirea lui k cu k+1 i se revine la pasul 3.

Metoda gradientului conjugat poate fi folosit i n rezolvarea sistemelor de ecuaii liniare .Mediul de calcul MATLAB

Pentru implementarea n Matlab a algoritmului metodei gradientului conjugat n versiunea Fletcher-Reeves este dezvoltat funcia fmingc utiliznd informaiile anterioare.

function [sol,gradi,nriter]=fmingc(f,gradf,f1,v,eps)

% Algoritmul de calcul al punctului de minim prin metoda gradientului conjugat in versiunea Fletcher-Reeves pentru o functie neliniara f(v) de n variabile

% f este functia de variabilele independente v, care reprezinta o functie definita de utilizator,

% gradf este gradientul functiei f, care reprezinta o functie definita de utilizator,

% f1 este functia de cautare liniara a pasului optimal, reprezentand de asemenea o functie definita de utilizator,

% v reprezinta o matrice coloana continand cele n valori initiale ale vectorului variabilei independente v,

% eps este precizia dorita.

% sol este solutia problemei de optimizare (valoarea optimala a lui v),

% gradi este gradientul functiei obiectiv, calculate in punctul v = sol,

% nriter este numarul de iteratii pana la obtinerea lui sol.

global p1 d1;

n=size(v);

nriter=0;

% Calculul gradientului initial:

df=feval(gradf,v);

% Bucla principala:

while norm(df)>eps,

nriter=nriter+1;

df=feval(gradf,v);

d1=-df;

% Bucla interna:

for intern=1:n,

p1=v;

% Precizia cautarii liniare = 0.00005; poate fi imbunatatita prin scadere

% Parametrii cautarii liniare:

options=optimset('Display','off','TolX',.00005);

pas=fminbnd(f1,-10,10,options);

% Calculul noii valori a lui v:

v1=v+pas*d1;

% Salvarea gradientului anterior:

aaa=df;

dfa=aaa';

% Calculul noului gradient:

aaa=feval(gradf,v1);

df=aaa';

% Actualizare v si d:

aaa=d1;

d=aaa';

v=v1;

% Formula specifica metodei Fletcher-Reeves:

betta=(df'*df)/(dfa'*dfa);

d1=-df+betta*d;

end

% Solutia:

sol=v1; % valoarea optimala a variabilei v

% Numarul de iteratii este in variabila nriter

gradi=df; % gradientul functiei obiectiv in punctul v optimal

disp('Numarul de iteratii='); disp(nriter);

disp('Solutia optimala='); disp(v1);

disp('Gradientul in solutia optimala='); disp(df);

Observaii:

1. Funcia fmincg necesit definirea de ctre utilizator a trei funcii Matlab.

2. n cadrul acestei funcii a fost utilizat funcia Matlab fminbnd pentru rezolvarea problemei de optimizare unidimensional avnd ca soluie valoarea optimal a pasului de cutare. Aceast funcie este nsoit de funcia Matlab optimset, prezentat, prin care sunt stabilite valorile unor parametri afereni funciei fminbnd.3. n cazul problemelor de optimizare la care funcia obiectiv este complicat trebuie modificat corespunztor funcia de cutare liniar a pasului

optimal.

Funcia fmincg va fi aplicat n cele ce urmeaz n cadrul exemplului urmtor:S se rezolve prin metoda gradientului conjugat n versiunea Fletcher-Reeves urmtoarea problem de optimizare:

Soluie: --Pentru nceput, este creat un fiier funcie, cu numele f91.m, pentru definirea funciei obiectiv:

**************************************************************************

function f=f91(v)

% Expresia functiei obiectiv:

f=0.5*(v(1)^4-16*v(1)^2+5*v(1))+0.5*(v(2)^4-16*v(2)^2+5*v(2));

**************************************************************************--Apoi, este creat fiierul funcie gradf91.m, destinat definirii gradientului funciei obiectiv:

**************************************************************************

function gradf=gradf91(v)

% Expresia gradientului functiei obiectiv:

gradf=zeros(size(v));

gradf(1)=0.5*(4*v(1)^3-32*v(1)+5);

gradf(2)=0.5*(4*v(2)^3-32*v(2)+5);

**************************************************************************--n continuare este creat i al treilea fiier funcie, numele fiierului fiind f191.m, pentru definirea funciei de cutare liniar utilizat n obinerea valorii optimale a pasului de cutare:

**************************************************************************

function pasn=f191(pas)

% Expresia functiei de cautare liniara utilizata in obtinerea valorii optimale a pasului de cautare:

global p1 d1;

q1=p1+pas*d1;

pasn=feval('f91',q1);

**************************************************************************

--n final, este executat urmtoarea secven de program Matlab care apeleaz funcia fmingc, rezultatele fiind prezentate imediat:

global p1 d1;

v0=[.6 .6];

x1=fmingc('f91','gradf91','f191',v0,.000005);

Rezultatul este afiat:

>> Numarul de iteratii=

2

Solutia optimala=

-2.9035 -2.9035

Gradientul in solutia optimala=

1.0e-006 *

-0.3815

-0.3815

Implementarea paralel al gradientului conjugatn cadrul fiecrei iteraii al algoritmului gradientului conjugat, se ia n calcul un singur vector-matrice. Acest calcul crete performana calculul gradientului conjugat i poate fi implementat n calcul paralel.

Pentru nceput matricea (A) i vectorul (x) se mpart de a lungul liniilor procesorilor multiplii ca n figura:

Exemplu de distribuie al A i x pe 4 procesoare

Algoritmul pentru calculul distribuit este: Pentru fiecare procesor j se distribuie x(j) i se calculeaz y(j)=A(j,:)*x

Dac avem la dispoziie un numr de n procesoare, este normal s lsm al i -lea procesor s controleze al i -lea component al vectorului n cauz, care este x(t), s(t) i g(t). Produsele acestor vectori sunt calculate prin permiterea procesorului al i -lea s calculeze produsul componentei al i -lea i adunarea sumelor pariale de-a lungul a 3 procesoare. Acesta este o acumulare de nod singular i necesit timp proporional cu diametrul reelei de interconectare. Dup aceasta valorile calculate sunt transmise la toate procesoarele. Presupunem acum c fiecare procesor a primit intrrile de la diferite rnduri ale lui A. Apoi vectorul matrice Ax(t) poate fi calculat prin transmiterea vectorului x(t) (transmisie multinod) i permiterea procesorului al i -lea s calculeze produsul lui x din rndul al i -lea al lui A. ntre timp, procesorul al i -lea poate calcula [A]jixi pentru fiecare j i aceste calcule pot fi propagate la fiecare procesor j.

Dac avem mai puin de n procesoare la dispoziie calculul este asemntor cu excepie c avem mai multe componente, mai multe rnduri pentru matricea A asignat pentru fiecare procesor.

Comunicarea colectiv

Operaiile colective implic un grup de procese. Pentru execuie, operaia colectiv trebuie apelat de toate procesele, cu argumente corespondente. Procesele implicate sunt definite de argumentul comunicator, care precizeaz totodat contextul operaiei.

Multe operaii colective, cum este difuzarea, presupun existena unui proces deosebit de celelalte, aflat la originea transmiterii sau recepiei. Un astfel de proces se numete rdcin. Anumite argumente sunt specifice rdcinii i sunt ignorate de celelalte procese, altele sunt comune tuturor.

Operaiile colective se mpart n mai multe categorii:

-Sincronizare de tip barier a tuturor proceselor grupului

-Comunicri colective, n care se include:

-Difuzarea

-Colectarea datelor de la membrii grupului la rdcin

-Distribuirea datelor de la rdcin spre membrii grupului

-Combinaii de colectare / distribuire (allgather i alltoall)

-Calcule colective:

-Operaii de reducere

-Combinaii reducere / distribuire

Calculul multithread n MatlabMulte din funciile Matlab permit calculele multithread i ne furnizeaz performane nbuntite la sistemele multiprocesor. Aceste funcii includ operaii de algebr liniar care apeleaz librria BLAS (de exemplu nmulirea matricilor). Exemplul urmtor de demonstreaz performanele nbuntite folosind calcul cu 2 threaduri.

Configurarea calculului multithread din meniul Preferences (Matlab R2007b):

Pentru performae optime este recomandat acceptarea valorii oferite ca valoare de pornire.Msurarea nbuntirii performanei pentru o operaie singular:

Acest exemplu folosete dou threaduri -definite n variabila numThreads pentru operaie de nmulire matrice cu matrice. Se poate continua experimentul prin nbuntirea i mai mult al performanei prin mrirea numrului threadurilor, dar acest lucru necesit mai multe procesoare.Primul pas este determinarea parametrilor i generarea de date aleatorii n variabilele A i B:

numThreads=2; % Numarul de threads pentru test

dataSize=500; % Marimea datelor de test

A=rand(dataSize,dataSize); % Matrice patratica randomB=rand(dataSize,dataSize); % Matrice patratica randomAl doilea pas este setarea numrului de threads la 1 i timpul:

oldstate = maxNumCompThreads(1);

C=A*B;

% setam numarul de threads si aflam timpul necesar pentru eltic;

C=A*B;

time1=toc;

fprintf('Time for 1 thread = %3.3f sec\n', time1);

Dup rulare obinem rezultatul:

Timpul obinut este de 0,091 secunde.Urmtorul pas este setarea numrului de threads la numThreads i timpul. Putem mri acest numr dac avem la dispoziie mai mult de 2 procesoare:

maxNumCompThreads(numThreads);

tic;

C=A*B;

timeN=toc;

fprintf('Time for %d threads = %3.3f sec\n', numThreads, timeN);

Dup cum se observ am obinut un timp de 0,068 secCalculm nbuntirea performanei:

Msurarea nbuntirii performanei pentru operaii multiplii:

Acest exemplu ilustreaz creterea performanei pentru funcii multiple. Folosete funcia ajuttoare runAndTimeOps pentru a calcula cteva rulri.

Se introduce n command line:

runAndTimeOpsdup care se creaz fiierul:

**************************************************************************

function [meanTime names] = runAndTimeOps% Masoara timpii necesari ptr un numar de operatii si returneaza timpii% necesari pentru ele si denumirea lor% Setarea parametrilornumRuns = 10; % Numarul de rulari mediidataSize = 500; % Marimea datelor de testx=rand(dataSize,dataSize); % Matrice aleatoare patratica% Inmultirea matricilor (*)func=1; % Initializarea functiei countertic;for i = 1:numRuns y=x*x; % Apelarea functieiendmeanTime(func)=toc/numRuns; % Imparte timpul cu numarul de rularinames{func}='inmult.'; % Memoreaza string care descrie functiafunc=func+1; % Incrementeaza functia counter% Impartirea matricilor (\)tic;for i = 1:numRuns y=x\x(:,1); % Apelarea functieiendmeanTime(func)=toc/numRuns; % Imparte timpul cu numarul de rularinames{func}='impart.'; % Memoreaza string care descrie functiafunc=func+1; % Incrementeaza functia counter% Functia sinus cu argumentul in radianitic;for i = 1:numRuns y=sin(x); % Apelarea functieiendmeanTime(func)=toc/numRuns; % Imparte timpul cu numarul de rularinames{func}='sinus'; % Memoreaza string care descrie functiafunc=func+1; % Incrementeaza functia counter% Ridicarea la putere a elementelor unui vector/matricetic;for i = 1:numRuns y=x.^x; % Apelarea functieiend

meanTime(func)=toc/numRuns; % Imparte timpul cu numarul de rularinames{func}='.^'; % Memoreaza string care descrie functiafunc=func+1; % Incrementeaza functia counter% radacina patraticafor i = 1:numRuns y=sqrt(x); % Apelarea functieiendmeanTime(func)=toc/numRuns; % Imparte timpul cu numarul de rularinames{func}='sqrt'; % Memoreaza string care descrie functiafunc=func+1; % Incrementeaza functia counter% Inmultirea element cu element (.*)tic;for i = 1:numRuns y=x.*x; % Apelarea functieiendmeanTime(func)=toc/numRuns; % Imparte timpul cu numarul de rularinames{func}='.*'; % Memoreaza string care descrie functiafunc=func+1; % Incrementeaza functia counter

**************************************************************************

Apelm funcia. Setm numrul de threads i timpul.maxNumCompThreads(1);

% Setarea numarului de threads[time1thread functionNames]=runAndTimeOps; % Timpii operatiunilor

Setm numrul de threads la numThreads i timpul din nou:maxNumCompThreads(numThreads);

% Setarea numarului de threads[timeNthreads functionNames]=runAndTimeOps; % Timpii operatiunilor

Restabilim numrul de threads avute nainte de rularea exemplului:maxNumCompThreads(oldstate);%reinitializeaza numarul de calcule cu valorile precedente

Calculm nbuntirile de performan:

speedup=time1thread./timeNthreads; % Viteza tuturor functiilor

Desenm pe un grafic nbuntirile:

bar(speedup); % Deseneaza ca bare vitezele pentru toate operatiile

title(['Cresterea Performantei cu ' int2str(numThreads) ' (nr. de threads)pe o arie de ' int2str(dataSize) 'x' int2str(dataSize)]);

ylabel('Cresterea Performantei');

set(gca, 'XTickLabel', functionNames);

ylim([0 2.25]); % valorile maxime ale axei Y

xlim([0 length(functionNames) + 1]); % valorile maxime ale axei X

grid;Rezultatul fiind:

Concluzie

De mai bine de 50 de ani, algoritmii de gradient conjugat au fost obiectul unor cercetri i analize teoretice i computaionale intensive. i n prezent, aceti algoritmi continu s reprezinte o component important a metodelor de optimizare fr restricii, dovedindu-i caracterul lor de inepuizabilitate. n aceast lucrare, am prezentat cteva probleme deschise, care constituie subiecte de meditaie i cercetare pentru proiectarea i implementarea n programe de calcul de noi algoritmi eficieni de gradient conjugat.

_1325362449.unknown

_1325431146.unknown

_1325362403.unknown


Recommended