% sr1t.m -- SR1-based Quasi-Newton script, Trust Region method
%
% Script applies SR1-based TR method to a given function, 
% using workspace variables described below:
%   fname  - string containing name of function 
%   dfname - string containing name of gradient 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
%
% Uses approximately exact solutions to the TR subproblem.

% (c) 2001, Philip D. Loewen
% 25 Feb 01  --  Original
% 12 Mar 01  --  Modified TR radius update rule

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

Delta = 0.001;       % 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(['SR1/Trust Region minimization of function "',fname,'", ',...
      'using gradient "',dfname,'":']);

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

% Count flops for added interest
fsf = flops;

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

% Initialize approx Hessian
Id = eye(max(size(x0)));
Hc = Id*abs(fc);    
% Build output line
disp(['  0', sprintf('  %12.4e',[xc',fc]),...
      sprintf('  %10.2e',Delta)]);

for k=1:Kmax,
  
  % Evaluate Newton point.
  pN = - Hc\gc';
  predredN = -(gc*pN + 0.5*pN'*Hc*pN); % Model reduction for Newton Pt

  
  lambda1 = min(eig(Hc));
  if pN'*pN <= Delta^2 & lambda1>-10*eps,
    xt = xc + pN;    % Convex model, Newton point in trust region, done!
  else
    % Newton iteration to find suitable lam>-lambda1
    lam = -lambda1 + 0.25*abs(lambda1) +100*eps;
    plam = -(Hc+lam*Id)\gc';
    xsnum = 0;
    HardCase = 0;
    while abs(plam'*plam-Delta^2)>0.01*Delta^2 & xsnum<5,
      xsnum = xsnum + 1;
      R = chol(Hc+lam*Id);
      qlam = R'\plam;
      np2 = plam'*plam;
      lam = lam + np2/(qlam'*qlam)*(sqrt(np2)-Delta)/Delta;
      if lam < -lambda1, HardCase=1; break; end
      plam = -(Hc+lam*Id)\gc';

    end
    % Check for Hard Case and use Cauchy Point then
    % Check if xsnum=5 and think about this.
    if 0,
      disp(' ');
      HCS = ''; if HardCase, HCS=' (Hard Case)'; end
      disp(['     xsnum=',int2str(xsnum),...
          ', lambda1=',num2str(lambda1),...
          ', lam=',num2str(lam),...
          ', |p(lam)|=',num2str(norm(plam)),...
          ', Delta=',num2str(Delta),...
          HCS,':']);
    end
    xt = xc + plam;
  end
  
  % Trust-Region subproblem has been "solved" above:
  % solution is called xt for "x-trial".
  % Find function and gradient values.
  ft = eval([fname,'(xt)']);
  gt = eval([dfname,'(xt)']); 
  
  dx = xt - xc;   % Change in x (a column vector).
  df = ft - fc;   % Change in f-values (should be negative).
  dg = gt - gc;   % Change in gradient values (a row vector).

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

  %  disp(' ');
  %  disp([sprintf('%3d  actred = %10.2e, predred = %10.2e, eig(Hc) = ',...
  %                                                 [k,actred,predred]),...
  %        sprintf(' %10.2e',eig(Hc)')] );
    
  % Update the current point
  if RR>=RRT,
    xn = xt;    % Good reduction, accept trial point.
    fn = ft;
  else
    xn = xc;    % Poor reduction, retain current point.
    fn = fc;
  end;
  gn = eval([dfname,'(xn)']);
  
  % Evaluate stopping criteria.
  e1 = max(abs(pN') ./ max([abs(xc');Zeps*ones(size(xc'))])); % Rel length of Newton step
  e2 = max( (abs(gn') .* abs(xn)) / max([abs(fn),Zeps]) );    % Max comp of grad(f)

  % Build and display step report
  
  disp([sprintf('%3d',k), sprintf('  %12.4e',[xn',fn]),...
        sprintf('  %10.2e',Delta),...
        sprintf('  %6.4f',norm(dx)/Delta),sprintf('  %7.3f',RR),...
        sprintf('  %8.1e',e1),...
        sprintf('  %8.1e',e2)]);

  xc = xn;                     % Update the x-vector
  fc = fn;                     % Update the function value
  gc = gn;                     % Update the gradient
  
  % Update the trust region radius:
  if RR<0.25, Delta = 0.25*norm(dx,2); end;
  if RR>0.75 & norm(dx,2)>((1-eps)*Delta),
    Delta = min([2*Delta;Dmax]);
  end;
  
  % Act on stopping criteria.
  stopcode = zeros(2,1);
  stopcode(1) = (e1<eta1 & predredN>0);
  stopcode(2) = (e2<eta2);
  if stopcode(1)|stopcode(2), break; end;
  
  % Use SR1 to update the approximate Hessian
  % Recall TR subprob solution produced trial step dx,
  % gradient difference dg.
  y_Hs = dg' - Hc*dx;
  if (y_Hs'*y_Hs>eps & (dx'*y_Hs)^2 >= 1.e-14*(dx'*dx)*(y_Hs'*y_Hs)),
    Hc = Hc + (y_Hs*y_Hs')/(dx'*y_Hs);
  else  
    disp('     Hessian update skipped.');
  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 stopcode(1), 
  ol = 'e1 < eta1'; 
  if stopcode(2),
    ol = [ol,' and '];
  end;
end;
if stopcode(2),
  ol = [ol,'e2 < eta2'];
end; 
if k==Kmax,
   ol = 'iteration limit reached';
end;
disp(['Stopped by:  ',ol,'.']);

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

