%===============================================
% Description : Steepest descent demonstration
%               Minimizes a quadratic function
%
%         J(x) = (1/2)x^T A x - b^T x
%
% Author : C. Alex Simpkins
% Adapted from code by Prof. Thomas R. Bewley, UCSD
%
%
%===============================================


%we want to solve the simultaneous equations for x1 and x2
%  4*x1 + 1*x2 = 3
%  1*x1 + 5*x2 = 3



% A must be positive definite...
% We want Ax=b form, so we build A from the above equation by inspection, 
% and b as well...

A = [4 1;...
     1 5];
	 
b = [3; 3];

%so we can plot our objective function, compute it and store the results...
[J,X,Y] = precompute_cost(A,b,100);

%clear extraneous values from any previous runs...
clear res_save x_save x; 


%plot the objective function as a surf...
figure(1);
pcolor(X,Y,J); shading interp; 
axis([-3 6 -3 6]);
xlabel('x_1');
ylabel('x_2');
hold on;

%the actual algorithm...
%-----------------------
%tolerance and maximum number of iterations...
epsilon = 1e-6; 
itmax = 1000;

%our initial guess...
x=zeros(size(b));
x=[-100; -100];


for iter=1:itmax
    r=(b-A*x);								 % determine gradient
    res=r'*r;								 % compute residual
    res_save(iter)=res;						 % save some stuff for later
    x_save(:,iter)=x; 
    if(res<epsilon | iter==itmax),break; end % exit yet?
    alpha=res/(r'*A*r);						 % compute alpha
    x=x+alpha*r;							 % update x
end
%-----------------------


disp('Number of iterations is'); iter

%just plot the parameters over the iterations...
for i=1:(length(x_save(1,:))-1),
    plot([x_save(1,i),x_save(1,i+1)], [x_save(2,i),x_save(2,i+1)],'y')
    %pause;
end

figure;
plot(res_save(1,:),'r')
xlabel('iteration number')
ylabel('error=(y-yp)^2')
title('Residuals vs. iteration number for gradient descent')