% cauchypt.m -- Trust Region/Cauchy Point script for minimization
%
% Script applies TR-based Steepest Descent method to a given function, 
% using workspace variables described below:
%   fname  - string containing name of function 
%   dfname - string containing name of gradient function
%   hfname - string containing name of Hessian function
%   x0     - column vector containing starting point
% Stopping test parameters are defined near the top, and can be
% changed without harming what remains.  These are
%   eta1   - relative error among components of the input vector x
%   eta2   - tolerance for the gradient being nearly zero
%   Kmax   - maximum number of iterations

% (c) 2001, Philip D. Loewen
% 23 Feb 01  --  Original
% 10 Mar 01  --  Modified Trust Region Radius Update

eta2 = 10000*eps;  % gradient-based termination tolerance.
Kmax = 500;          % Students were asked to use 25
RRT  = 1/5;        % Reduction Ratio Threshhold ("eta" in notes).

Delta = 0.1;       % Initial value for trust region radius.
Dmax = 100;        % Absolute maximum trust region radius.

Zeps = eps^(3/4);  % Prevents division by 0 in stopping tests.

% Set up table titles
disp([' ']);
disp(['Trust Region/Cauchy Point minimization of function "',fname,'",']);
disp(['using gradient "',dfname,'", Hessian "',hfname,'":']);

xstring = '';
for jj=1:length(x0),
  xstring = [xstring,'x_k(',int2str(jj),')        '];
end;
disp([' ']);
disp(['   k     ',xstring,'f(x_k)        Delta   |dx|/Delta   RR']);

% Count flops for added interest
fsf = flops;

% Set current point, function, and gradient.
xc = x0;
fc = eval([ fname,'(xc)']);
gc = eval([dfname,'(xc)']);

% Build output line
disp(['   0', sprintf('  %12.4e',[xc',fc]),...
      sprintf('  %10.2e',Delta)]);

for k=1:Kmax,

  % Evaluate current Hessian, test concavity along gradient direction.
  Hc = eval([hfname,'(xc)']);
  gHg = gc*Hc*gc';
  
  if gHg>5*eps,
    % 1D function is convex enough:  
    % Find Unconstrained Steepest Descent Minimizer for Model Function.
    pSD = -(gc*gc')/gHg*gc';
    pSDlen = (pSD'*pSD)^0.5;
    if pSDlen <= Delta,
      xCP = xc + pSD;
    else
      xCP = xc + Delta*pSD/pSDlen;
    end;
  else
    % 1D slice of model function in gradient direction is concave:
    % SD minimizer lies on TR boundary, in direction of -g'
    xCP = xc - Delta*gc'/(gc*gc')^0.5;
  end;
  
  % Trust-Region subproblem has been "solved" by Cauchy point.
  % Find function and gradient values.
  fCP = eval([fname,'(xCP)']);
  
  dx = xCP - xc;   % Record change in x-vector.
  df = fCP - fc;   % Record change in f-values (should be negative).

  actred = -df;   % Reduction observed on objective function.
  predred= -gc*dx - 0.5*dx'*Hc*dx;  % Reduction predicted by model.
  RR     = actred/predred;   % Reduction Ratio.

  % Move to new point only if Reduction Ratio is large
  if RR>RRT,
    xc = xCP;
    fc = fCP;
    gc = eval([dfname,'(xc)']); 
  end;
 
  % Evaluate gradient-based stopping criterion.
  e2 = max( (abs(gc') .* abs(xc)) / max([abs(fc),Zeps]) );  % Max rel chg in df

  % Build output line
  disp([sprintf('%4d',k), sprintf('  %12.4e',[xc',fc]),...
        sprintf('  %10.2e',Delta),...
        sprintf('  %7.5f',norm(dx)/Delta),sprintf('  %8.5f',RR)]);
  
  % Update the trust region radius:
  if RR<0.25, Delta = 0.25*norm(dx); end;
  if RR>0.75 & norm(dx,2)>((1-eps)*Delta),
    Delta = min([2*Delta;Dmax]);
  end;
  
  % Act on stopping criteria.
  if e2<eta2, break, end
  
end;    % Finished main iteration.

totflops = flops-fsf;

% Print results.
disp([' ']);
disp(['Best point:  x'' =',sprintf(' %22.16e',xc),'.']);
disp(['Best value:  ',fname,'(x) =',sprintf(' %22.16e',fc),'.']);
disp(['Step count:  ',int2str(k),' Trust-Region passes.']);

ol = '';
if (e2<eta2), 
  ol = [ol, 'e2 < eta2']; 
end;
if isempty(ol),
  ol = 'iteration limit reached';
end;
disp(['Stopped by:  ',ol,'.']);

disp(['Flop count:  ',int2str(totflops),' total--average ',int2str(totflops/k),' per step.']);
disp(' ');

