% minimize f(x)=.5(alpha*x_1^2+x_2^2)
%% Run this for different alpha values and different starting points
hold off
clear all
clf
alpha=16;   % run for alpha values  2,4,8,16, etc...
Q=[alpha 0;0 1];   % Hessian
eigQ=eig(Q);
ratio2=4*prod(eigQ)/(sum(eigQ))^2;
K=(alpha-1)/(alpha+1);
disp(['decrease constant K is:  ',num2str(K)]);
X= [-11: .4: 11];
Y= [-11: .4: 11];
v=  alpha*X.^2 + Y.^2;
[X,Y] = meshgrid([X]);
Z = .5*(alpha*X.^2 + Y.^2);
w=[v(1) v(round(length(v)/4) : end-round(length(v)/4) )   v(end)];
[C,H]=contour(X,Y,Z,w);
hold on
xnew=[10/alpha 10]';    % starting points - worst case convergence
xs=[];
%xnew=[10 10]';    % starting point
%%xnew=[10 0]';    % starting point  - best case
gradnew=Q*xnew;
qnew=.5*xnew'*Q*xnew;
disp(['iterate ',num2str(1),' is:   ',num2str(xnew') '  with value: ',...
               num2str(qnew) ]);
plot(xnew(1),xnew(2),'xk')
disp(['initial point: ',num2str(xnew')])
%disp(['hit any key to continue'])
%pause
for i=1:100,
%% update
x=xnew;
xs=[xs x];
q=qnew;
grad=gradnew;
tstep=norm(grad)^2/(grad'*Q*grad);
xnew=x-tstep*grad;
gradnew=Q*xnew;
disp(['cos angle between grads is:  '...
      ,num2str( gradnew'*grad/(norm(gradnew)*norm(grad)) )]);
%plot(xnew(1),xnew(2),'xk')
plot(xs(1,:),xs(2,:),xnew(1),xnew(2),'xk');
qnew=.5*xnew'*Q*xnew;
disp(['iterate ',num2str(i+1),' is:   ',num2str([xnew'])...
      '  with value: ', num2str(qnew) ]);
ratio1=(grad'*grad)^2 / ( grad'*Q*grad*grad'*(Q\grad) ) ;
disp(['exact guaranteed decrease is: ',num2str(...
       q- q*ratio1)]);
disp(['guaranteed decrease using K is: ',num2str(K^2*q)])
disp(['enter to continue'])
keyboard
%pause(.01)     % change to pause to view decrease
end
hold off

