% quadratic case - steepest descent
clear
X= [1 1; 0 1];Y=[-1; 1];
BOPT=inv(X)*Y
for i=1:41
   xx(i)=-4+(i-1)*0.1;
   for j=1:41
      yy(j)=-1+(j-1)*0.1;
      BB=[xx(i);yy(j)];
      Z(i,j) = (Y-X*BB)'*(Y-X*BB);
   end
end
figure(1); clf;
contour(xx,yy,Z,20)
xlabel('B_1')
ylabel('B_2')
hold on
clear xx; clear yy   
xtickformat('%0.1f')
ytickformat('%0.1f')

B=[0;0];
for k=1:10
   yy(k)=B(2);xx(k)=B(1);
   gr=-X'*(Y-X*B);
   norm1=gr'*gr; norm2=(X*gr)'*(X*gr);
   B=B-norm1/norm2*gr
   s(k)=0.5*(Y-X*B)'*(Y-X*B)
end
plot(xx,yy,'lineWidth',2)
figure(2); clf;
semilogy(s, 'linewidth', 1.0)
xlabel('iteration');
ylabel('S=(Y-XB)^T(Y-XB)')
ytickformat('%0.2f')
grid on

% 
% fig=figure(2);
% print(fig,fullfile(pwd, 'Figure6_10.tif'),'-dtiff','-r600')
% 
% fig=figure(1)
% print(fig,fullfile(pwd, 'Figure6_11.tif'),'-dtiff','-r600')
%