0001
0002
0003
0004
0005
0006
0007
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
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
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
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
0120 end
0121
0122
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