function BFGS

tol=10^(-5);

rho =1/8;
c=10^(-4);

x=[1.2; 1.2];
[y,gy,hy] =rosenbrock(x);

Xplot=x;
counter=0;

B=eye(2);
while norm(gy)>tol
    
    p=-B\gy;
    gyold=gy;
    
    alphak=linesearch(@rosenbrock,p,x,rho,c);
    x=x+alphak*p;
    [y,gy,hy] = rosenbrock(x);
   
    
    s= alphak*p;
    y= gy-gyold;
    
    B=B-(B*s*s'*B)/(s'*B*s)+(y*y')/(y'*s);
    
    
    Xplot=[Xplot x];
    counter=counter+1;
end

counter

close all
[X,Y]=meshgrid(-2:0.01:2, -1:0.01:3);

Z=100*(Y-X.^2).^2+(1-X).^2;

v=logspace(0.5,3,30);
contour(X,Y,Z,v);
hold on

for j=2:counter
    plot([Xplot(1,j),Xplot(1,j-1)], [Xplot(2,j), Xplot(2,j-1)],'-o')
end
hold off

end

%%%%%%%%%%%%%%%%%%%%%%%%


function alphak = linesearch(f,d,x,rho,c)

% Backtracking line search
% See Algorithm 3.1 on page 37 of Nocedal and Wright
% Input Parameters :
% f: MATLAB file that returns function value
% d: The search direction
% x: previous iterate
% rho :- The backtrack step between (0,1) usually 1/2
% c: parameter between 0 and 1 , usually 10^{-4}
% Output :
% alphak: step length calculated by algorithm

alphak = 1;
[fk, gk] = feval(f,x);
xx = x;
x = x + alphak*d;
fk1 = feval(f,x);
while fk1 > fk + c*alphak*(gk'*d)
  alphak = alphak*rho;
  x = xx + alphak*d;
  fk1 = feval(f,x);
end

end

