newton.m

求解无约束最小化问题的牛顿方法 matlab程序

% newton.m -- Pure Newton script for unconstrained minimization
%
% Script applies Newton's 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
%   eta3   - tolerance for the improvement in function value
%   Kmax   - maximum number of iterations

% (c) 2001, Philip D. Loewen
% 21 Feb 01  --  Original
% 23 Feb 01  --  Cosmetics
% 24 Feb 01  --  Cosmetics
% 11 Mar 01  --  Big overhaul, new stopping tests, typical values

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

vsn  = sqrt(realmin);       % Very Small Number
vsv  = vsn*ones(size(x0));  % Very Small Vector 
vbv  = ones(size(x0))./vsn; % Very Big Vector

% Set up table titles
disp([' ']);
disp(['Pure Newton 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)      dx(rel)   grad(rel)',...
      ' pred obj dec']);

% Count flops for added interest
fsf = flops;

% Set current point, old point, old function value.
xc = x0;
fc = eval([ fname,'(xc)']);
xo = 2*(abs(x0) + ones(size(x0)));
fo = 2*abs(fc);

% Main Loop of Newton's Method:
for k=0:Kmax-1,
  fc = eval([ fname,'(xc)']);
  gc = eval([dfname,'(xc)']);
  Hc = eval([hfname,'(xc)']);
  
  % Evaluate stopping criteria.
  xtyp = (abs(xc) + abs(xo) + vsv)/2;
  ftyp = (abs(fc) + abs(fo) + vsn)/2;

  e1 = max(abs(xc-xo) ./ xtyp);                % Rel chg in x
  e2 = max( (abs(gc') .* abs(xtyp)) / ftyp );  % Max rel chg in grad(f)
  e2 = max([e2, norm(gc)*norm(xtyp)/ftyp]);    % Alt form of above
  e3 = 0.5*gc*Hc*gc'/ftyp;                     % Rel chg in f pred by Newton
  
  % Build and print output line
  ol = [sprintf('%3d',k), ...  
        sprintf('  %12.4e',xc'), ... 
        sprintf('  %12.4e',fc)];
  if fc-fo>10*eps,
    ol = [ol,'!'];
  else
    ol = [ol,' '];
  end;   % "!" = "objective increased!"
  ol = [ol,sprintf(' %9.1e  %9.1e  %9.1e',[e1,e2,e3])];
  disp(ol);
  
  tests = [(e1<eta1), (e2<eta2), abs(e3)<eta3];
  if sum(tests), break, end

  xo = xc;                     % Save old x-vector for stopping tests.
  fo = fc;                     % Save old f-value  for stopping tests.
  
  p  = - Hc \ gc';             % Full Newton step from current point xc.
  xc = xc + p;                 % Update.
  
end;    % Finished main iteration.

totflops = flops-fsf;

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

ol = '';
while sum(tests),
  tn = min(find(tests));
  if ~isempty(ol),
    ol = [ol,', '];
  end
  ol = [ol,'test ',num2str(tn)];
  tests(tn)=0;
end;
if isempty(ol),
  ol = 'iteration limit reached';
end;
disp(['Stopped by:  ',ol,'.']);

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