nad.m

一个用拟牛顿算法求解优化问题的程序

% nad.m -- Newton/Armijo/Descent script for unconstrained minimization
%
% Script applies lightly-modified Newton 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
% 21 Feb 01  --  Original
% 23 Feb 01  --  Cosmetics
% 24 Feb 01  --  Limited Armijo over-relaxation steps to 10 in SD case

eta1 = 100*sqrt(eps);
eta2 = 10000*eps;
Kmax = 100;       % Students were asked to use 25

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

% Armijo line search parameters
alpha = 1.0e-3;  % Slope multiplier for "sufficient decrease" condition
beta  = 0.9;     % Geometric factor for backtracking iteration

% Set up table titles
disp([' ']);
disp(['Newton/Armijo/Descent minimization of function "',fname,'",']);
disp(['using gradient "',dfname,'" and Hessian "',hfname,'":']);

xstring = '';
for jj=1:length(x0),
  xstring = [xstring,'x_k(',int2str(jj),')        '];
end;
disp([' ']);
disp(['  k     ',xstring,'f(x_k)        e_1           e_2         lambda']);

% Count flops for added interest
fsf = flops;

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

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

for k=1:Kmax,

  % Set line search direction.  First, find Newton step:
  s  = - Hc\gc';

  % Trade Newton step for Steepest Descent step in bad case:
  SD = 0;     % Steepest descent flag.
  m  = gc*s;  % Slope in direction s from xc; expect m<0.
  if (m > -eps^0.25*norm(gc,2)*norm(s,2)),
    % Newton step is rather poor descent direction,
    % so just follow steepest descent direction instead.
    s = -gc';
    m = gc*s;
    SD = 1;
  end;

  % Take an Armijo step in direction s.
  lam = 1;         % Stretch factor lambda
  while 1,         % Repeat forever:
    xn = xc + lam*s;                                 % New point
    fn = eval([ fname,'(xn)']);                      % Function value
    if ( (fn-fc) <= (alpha*lam*m)+eps ), break, end; % Test acceptability
    lam = lam*beta;                                  % Shrink lambda, retry
  end;
  
  % In SD case, try up to 10 expansion steps:
  if SD & (lam==1),
    for dummy=1:10,
      lam = lam/beta;               % Stretch step length a little.
      xt = xc + lam*s;              % Set trial point in direction s.
      ft = eval([ fname,'(xt)']);   % Function value at trial point.
      if ft>fn, break; end;         % If higher, quit looking (keep prev xn,fn).
      xn = xt;                      % Otherwise, use new xn,fn.
      fn = ft;
    end;
  end;

  % Get here with new point xn and fcn value fn in good shape.
  % Work out new gradient too.
  gn = eval([dfname,'(xn)']); 
  dx = xn - xc;   % Record change in x-vector.
  dg = gn - gc;   % Record change in gradient.

  % Evaluate stopping criteria.
  e1 = max(abs(dx') ./ max([abs(xc');Zeps*ones(size(xc'))])); % Rel chg in x
  e2 = max( (abs(gn') .* abs(xn)) / max([abs(fn),Zeps]) );  % Max rel chg in df

  % Build output line
  ol = sprintf('%3d',k);
  ol = [ol, sprintf('  %12.4e',xn')];
  ol = [ol, sprintf('  %12.4e',[fn,e1,e2])];
  ol = [ol, sprintf('  %8.6f',lam)];
  if SD, ol = [ol,' (sd)']; end;   % "sd" = "steepest descent"
  disp(ol);

  xc = xn;                     % Update the x-vector
  fc = fn;                     % Update the function value
  gc = gn;                     % Update the gradient
  Hc = eval([hfname,'(xc)']);  % Update the Hessian

  if (e1<eta1)|(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),' modified Newton steps.']);

ol = '';
if e1<eta1, 
  ol = 'e1 < eta1'; 
  if e2<eta2,
    ol = [ol,' and '];
  end;
end;
if e2<eta2,
  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(' ');