function [Results, Count]=GradientDescent(x0,M)
%[Results, Count]=GradientDescent(x0,M)
%   Gradient descent method f(x)=1-exp(-x'*M*x) with initial guess
%   x0. We return a table Results, where each line contains the
%   iteration count, the current minimizer guess, the function
%   value at this guess, and the 2-norm of the gradient at this
%   guess. Count is the number of iterations to reach convergence.
% 
%   ($RCSfile: GradientDescent.m,v $ $Revision: 1.1.1.1 $)

% Generate grid.
a=1;x=[-a:0.05:a];y=[a:-0.05:-a];
% Compute function values.
[YY,XX]=meshgrid(x',y');
for i=1:length(x)
  for j=1:length(y)
    ZZ(i,j)=MyFun(M,[XX(i,j);YY(i,j)]);
  end
end
% Plot function surface.
figure;surf(XX,YY,ZZ);shading interp; colormap(autumn)
xlabel('x');ylabel('y');view(2); axis('square')
hold on

offset=100;               % to plot lines on top of surface
fx0=MyFun(M,x0);          % function value in current value x0
gradx0=MyGrad(M,x0);      % gradient value in current value x0
absgrad=norm(gradx0);     % 2-norm of current gradient
Count=0;                  % iteration counter
Results=[Count x0' fx0, absgrad];  % store statistics
mu=1;                     % initialize leran step size
while absgrad>1e-3        % loop until convergence
Count=Count+1;            % increase step counter
  xneu=x0-mu*gradx0;      % compute new guess
  fneu=MyFun(M,xneu);     % and function value at new guess
  while fneu > fx0        % do not allow increase in function value
    mu=mu/3;              % decrease mu
    xneu=x0-mu*gradx0;    % compute next new guess
    fneu=MyFun(M,xneu);   % and corresponding function value
  end                     % new guess is accepted
  linie=[ x0 xneu];       % generate line coordinates for plot
  z=[fx0 fneu]+offset ;   % generate height values for plot
  line(linie(1,:), linie(2,:), z,'Marker', '+'); % plot line
  x0=xneu;                % set x0 to new guess
  fx0=fneu;               % update function value
  gradx0=MyGrad(M,x0);    % compute new gradient
  absgrad=norm(gradx0);   % and its 2-norm
  Results=[Results;[Count x0' fx0, absgrad]]; % store statistics
  pause(0.005)            % wait a bit
  mu=2*mu;                % double learn step size
end                       % while loop completed

return

function f=MyFun(M,x)
  f=1-exp(-x'*M*x);
return 

function g=MyGrad(M,x)
  g=exp(-x'*M*x)*(M+M')*x;
return
