% Preconditioned Conjugate Gradient Method to Find the Minimum of the
% Rosenbrock Functions
%
% xstart :starting point.
% xact : current solution at one iteration.
% xnext : new solution after updating.
% step: recording of the moving path of the solution.
% 
% By Wei Dai Oct 30, 2011
clear all;
close all;
clc;


xstart=[-1 1]';
HH=zeros(2,2);
iteration=0;
delta=0.000000001;
xact=xstart;
step=zeros(2,5000);
step(1,1)=xstart(1);
step(2,1)=xstart(2);
dx1=zeros(2,1);
dx2=dx1;
g2=dx1;
g1=dx1;

%control parameters
pre=0; %use preconditioner or not
num=0; %use numerical line search or not
 
for i=1:5000
    iteration=iteration+1;
    x1=xact(1);
    x2=xact(2);
    f0 =  100*(x2-x1^2)^2+(1-x1)^2;
   dfx_dx1 = -400*x1*(x2-x1^2)-2*(1-x1);
   dfx_dx2 = 200*(x2-x1^2);
   d2fx_dx12 =1200*x1^2-400*x2+2;
   d2fx_dx1dx2 = -400*x1;
   d2fx_dx22 = 200;
   H=[d2fx_dx12 d2fx_dx1dx2;d2fx_dx1dx2 d2fx_dx22];
   if(pre==1)
      HH=[1/H(1,1) 0;0 1/H(2,2)];
   else
      HH=eye(2);
   end
   g2=[dfx_dx1 dfx_dx2]';
   if(iteration==1||mod(iteration,2)==0)
        beta = 0;
   else
        beta=g2'*HH*g2/(g1'*HH*g1);
   end
   
    
   dx2=HH*g2+beta*dx1;
  % alpha=g2'*dx2/(dx2'*H*dx2);
 alpha=g2'*HH*g2/(dx2'*H*dx2);
   
   if(num==1)
       xnext=xact-0.5*alpha*dx2;
       x1=xnext(1);
       x2=xnext(2);
       f1 =  100*(x2-x1^2)^2+(1-x1)^2;
       xnext=xact-alpha*dx2;
       x1=xnext(1);
       x2=xnext(2);
       f2 =  100*(x2-x1^2)^2+(1-x1)^2;   
   
       gama=(0.5^2*(f0-f2)+(f1-f0))/(-f0+2*f1-f2);
       xnext=xact-gama*alpha*dx2;
       x1=xnext(1);
       x2=xnext(2);
       fx =  100*(x2-x1^2)^2+(1-x1)^2;
       fprintf('%d %d %d %e %e\n',iteration,xnext,gama,fx);
   else
       xnext=xact-alpha*dx2;
       x1=xnext(1);
       x2=xnext(2);
       fx =  100*(x2-x1^2)^2+(1-x1)^2;
       fprintf('%d %d %d %e\n',iteration,xnext,fx);
   end
   
   step(1,i+1)=xnext(1);
   step(2,i+1)=xnext(2);
   if(norm(fx)<0.0000001)
        break;
   end
   xact=xnext;
   dx1=dx2;
   g1=g2;
end

%plot contour lines of f(x);
xx1=[-1.2:0.01:1.2];
xx2=[-0.3:0.01:1.2];
[XX1,XX2]=meshgrid(xx1,xx2);
f =  100.*(XX2-XX1.^2).^2+(1-XX1).^2;
figure;
contour(XX1,XX2,f,50);

axis image;
hold on;

plot(step(1,1:iteration+1),step(2,1:iteration+1),'-rs','LineWidth',2);
hold off;
title('Preconditioned Conjugate Gradient Method to Find the Minimum of the Rosenbrock Function','Fontsize',10);
