% 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 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(' ');
