Home > dal > newton.m

newton

PURPOSE ^

newton - a simple implementation of the Newton method

SYNOPSIS ^

function [xx,fval,gg,status]=newton(fun, xx, ll, uu, Ac, bc, tol,finddir, info, verbose, varargin)

DESCRIPTION ^

 newton - a simple implementation of the Newton method

 Syntax:
  [xx,fval,gg,status]=newton(fun, xx, ll, uu, Ac, bc, tol, finddir, info, verbose, varargin);

 Copyright(c) 2009 Ryota Tomioka
 This software is distributed under the MIT license. See license.txt

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 % newton - a simple implementation of the Newton method
0002 %
0003 % Syntax:
0004 %  [xx,fval,gg,status]=newton(fun, xx, ll, uu, Ac, bc, tol, finddir, info, verbose, varargin);
0005 %
0006 % Copyright(c) 2009 Ryota Tomioka
0007 % This software is distributed under the MIT license. See license.txt
0008 %
0009 function [xx,fval,gg,status]=newton(fun, xx, ll, uu, Ac, bc, tol, ...
0010                                     finddir, info, verbose, varargin)
0011 
0012 if isempty(verbose)
0013   verbose=0;
0014 end
0015 
0016 if iscell(fun)
0017   fcnHessMult = fun{2};
0018   fun = fun{1};
0019 else
0020   fcnHessMult = [];
0021 end
0022 
0023 if isempty(finddir)
0024   if isempty(fcnHessMult)
0025     finddir = @dir_chol;
0026   else
0027     finddir = @dir_pcg;
0028   end
0029 end
0030 
0031 
0032 n = size(xx,1);
0033 if verbose
0034   fprintf('n=%d tol=%g\n', n, tol);
0035 end
0036 
0037 
0038 cc = 0;
0039 step = nan;
0040 while 1
0041   [fval,gg,H,info]=feval(fun, xx, info, varargin{:});
0042   
0043   if verbose
0044     fprintf('[%d] fval=%g norm(gg)=%g step=%g\n',cc, fval, norm(gg),step);
0045   end
0046 
0047   if info.ginfo<tol % || norm(gg)<1e-3
0048     if verbose
0049       fprintf('Optimization success! ginfo=%g\n',info.ginfo);
0050     end
0051     ret = 0;
0052     status=archive('ret','cc','info');
0053     break;
0054   end
0055 
0056   if isstruct(H)
0057     H.fcnHessMult=fcnHessMult;
0058   end
0059   
0060   dd = feval(finddir, gg, H);
0061 
0062   [step,fval, info]=linesearch(fun, fval, xx, ll, uu, Ac, bc, dd, info, varargin{:});
0063 
0064   if step==0
0065     ret = -2;
0066     fprintf('[newton] max linesearch reached. ginfo=%g\n',info.ginfo);
0067     status=archive('ret','cc','info');
0068     break;
0069   end
0070 
0071   xx = xx + step*dd;
0072   
0073   cc = cc+1;
0074   if cc>1000
0075     ret = -1;
0076     fprintf('[newton] #iterations>1000.\n');
0077     status=archive('ret','cc','info');
0078     break;
0079   end
0080 end
0081 
0082 
0083 function dd = dir_chol(gg, H)
0084 R = chol(H);
0085 dd = -R\(R'\gg);
0086 
0087 %dd = pcg(H, -gg, max(1e-6,tol*0.01));
0088 
0089 function dd = dir_pcg(gg, Hinfo)
0090 [dd,dum1,dum2]=pcg(Hinfo.fcnHessMult{1}, -gg, 1e-2,...
0091                    length(gg),Hinfo.prec,[],[],Hinfo.fcnHessMult{2:end}, Hinfo);
0092 % [dd,dum1,dum2] = pcg(@(xx)fcnHessMult(xx,Hinfo), -gg, max(1e-6,tol*0.01), length(xx),Hinfo.prec);
0093 
0094 function [step,fval, info] = linesearch(fun, fval0, xx, ll, uu, Ac, bc, dd, info, varargin)
0095 
0096 Ip=find(dd>0);
0097 In=find(dd<0);
0098 step=min([1.0, 0.999*min((xx(In)-ll(In))./(-dd(In))), 0.999*min((uu(Ip)-xx(Ip))./dd(Ip))]);
0099 
0100 xx0    = xx;
0101 
0102 cc = 1;
0103 while 1
0104   xx = xx0 + step*dd;
0105   
0106   if ~isempty(Ac)
0107     bineq = all(Ac*xx<=bc);
0108   else
0109     bineq = true;
0110   end
0111   
0112   if bineq && all(ll<=xx) && all(xx<=uu)
0113     [fval,info] = feval(fun, xx, info, varargin{:});
0114 
0115     if fval<fval0
0116       break;
0117     end
0118   else
0119     % keyboard;
0120   end
0121   
0122   %fprintf('step=%g fval=%g (fval0=%g)\n',step, fval, fval0);
0123   
0124   step = step/2;
0125   cc = cc+1;
0126   if cc>30
0127     fval=fval0;
0128     step = 0;
0129     break;
0130   end
0131 end

Generated on Sat 22-Aug-2009 22:15:36 by m2html © 2003