% file name:  steepdesc_newt.m 
% This Matlab code implements Cauchy's steepest descent method 
% using Armijo stepsize rule.  
% It terminates when the norm of the gradient is below 10^(-6).
% The function value is read from the file "func.m".
% The gradient value is read from the file "grad.m".

% Read in inputs
  n=input('enter the number of variables n ');
  x=input('enter the initial column vector x ');

% Armijo stepsize rule parameters
  sigma = .1;
  beta = .5;
  obj=func(x);
  g=grad(x);
  k=0;                                  % k = # iterations
  nf=1;                                 % nf = # function eval.
  t1=cputime;
  ng=1;					% nf = # gradient eval.
  nh=0;

% Begin method
  while  norm(g) > 1e-6
    H = hessian(x);
    nh=nh+1;
    if min(eig(H)) > 0 
      d = -H\g;                 % Newton direction
    else
      d = -g;                   % steepest descent direction
    end

% convergence safeguard 
    if g'*d>-1e-6*g'*g | norm(d)>1e6*norm(g)
      d = -g
    end

    a = 1;	
    newobj = func(x + a*d);
    nf=nf+1;
    dirderiv = g'*d;			
    while newobj-obj > sigma*dirderiv*a
      a = a*beta;
      newobj = func(x + a*d);
      nf=nf+1;
    end
    if (mod(k,10)==1) fprintf('%5.0f %5.0f %12.5e %g %g \n',k,nf,obj,a,norm(g)); end
    x = x + a*d;
    obj=newobj;
    g=grad(x); 
    ng=ng+1;
    k = k + 1;

  end
  
  fprintf(' iter= %g  #func= %g  #grad= %g  #Hess= %g  cpu= %g \n',k,nf,ng,nh,cputime-t1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%file name:  func.m
%This routine evaluates the least-square trignometric function (ref. Spedicato).
%	f(x) = sum_{i=1}^n  g_i(x)^2
% g_i(x) = n - sum_{j=1}^n cos(x_j) + i(1-cos(x_i)) - sin(x_i)

function y = func(x)
n=size(x,1);    %n is the number of rows of the column vector x.

b = n - ones(1,n)*cos(x);
g = b*ones(n,1) + [1:n]' - [1:n]'.*cos(x) - sin(x);
y = norm(g)^2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%file name:  grad.m
%This routine evaluates the gradient of the least-square trignometric function 
%(ref. Spedicato).
%	f(x) = sum_{i=1}^n  g_i(x)^2
% g_i(x) = n - sum_{j=1}^n cos(x_j) + i(1-cos(x_i)) - sin(x_i)

function y = grad(x)
n=size(x,1);    %n is the number of rows of the column vector x.

b = n - ones(1,n)*cos(x);
g = b*ones(n,1) + [1:n]' - [1:n]'.*cos(x) - sin(x);
y = 2*(sum(g)*sin(x) + g.*([1:n]'.*sin(x) - cos(x)));