0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026 function [xx, status] = lbfgs(fun, xx, ll, uu, Ac, bc, info, opt, varargin)
0027
0028 opt=set_defaults(opt, 'm', 6,...
0029 'ftol', 1e-5, ...
0030 'maxiter', 0,...
0031 'max_linesearch', 50,...
0032 'display', 0,...
0033 'epsg', 1e-5,...
0034 'epsginfo', 1);
0035
0036 if ischar(fun)
0037 fun = {fun};
0038 end
0039
0040 nn = size(xx,1);
0041
0042 t0 = cputime;
0043
0044
0045 lm = repmat(struct('s',zeros(nn,1),'y',zeros(nn,1),'ys',0,'alpha',0),[1, opt.m]);
0046
0047 [fval,gg,info]=feval(fun, xx, info, varargin{:});
0048
0049
0050
0051 dd = -gg;
0052
0053 kk = 1;
0054 stp = 1/norm(dd);
0055
0056 bResetLBFGS = 0;
0057 ixend = 1;
0058 bound = 0;
0059 while 1
0060 fp = fval;
0061 xxp = xx;
0062 ggp = gg;
0063
0064
0065 [ret, xx,fval,gg,info,stp]=...
0066 linesearch_backtracking(fun, xx, ll, uu, Ac, bc, fval, gg, dd, stp, info, opt, varargin{:});
0067
0068 if ret<0
0069 fprintf('ginfo=%g\n',info.ginfo);
0070 break;
0071 end
0072
0073
0074 gnorm = norm(gg);
0075 if opt.display>1
0076 fprintf('[%d] xx=[%g %g...] fval=%g gnorm=%g step=%g\n',kk,xx(1),xx(2),fval,gnorm,stp);
0077 end
0078
0079 if info.ginfo<opt.epsginfo
0080 if opt.display>1
0081 fprintf('Optimization success! ginfo=%g\n',info.ginfo);
0082 end
0083 ret=0;
0084 break;
0085 end
0086
0087 if kk==opt.maxiter
0088 if opt.display>0
0089 fprintf('Maximum #iterations=%d reached.\n', kk);
0090 end
0091 ret = -3;
0092 break;
0093 end
0094
0095
0096 if opt.m>0
0097 lm(ixend).s = xx-xxp;
0098 lm(ixend).y = gg-ggp;
0099 ys = lm(ixend).y'*lm(ixend).s; yy = sum(lm(ixend).y.^2);
0100 lm(ixend).ys = ys;
0101 else
0102 ys = 1; yy = 1;
0103 end
0104
0105 bound = min(bound+1, opt.m);
0106 ixend = (opt.m>0)*(mod(ixend, opt.m)+1);
0107
0108
0109 dd = -gg;
0110
0111 jj = ixend;
0112 for ii=1:bound
0113 jj = mod(jj + opt.m -2, opt.m)+1;
0114 lm(jj).alpha = lm(jj).s'*dd/lm(jj).ys;
0115 dd = dd -lm(jj).alpha*lm(jj).y;
0116 end
0117
0118 dd = dd *(ys/yy);
0119
0120 for ii=1:bound
0121 beta = lm(jj).y'*dd/lm(jj).ys;
0122 dd = dd + (lm(jj).alpha-beta)*lm(jj).s;
0123 jj = mod(jj,opt.m)+1;
0124 end
0125
0126 stp = 1.0;
0127
0128 kk = kk + 1;
0129 end
0130
0131 status=struct('ret', ret,...
0132 'kk', kk,...
0133 'fval', fval,...
0134 'gg', gg,...
0135 'time', cputime-t0,...
0136 'info', info,...
0137 'opt', opt);
0138
0139
0140 function [ret, xx, fval, gg, info, step]...
0141 =linesearch_backtracking(fun, xx, ll, uu, Ac, bc, fval, gg, dd, step, info, opt, varargin)
0142
0143 floss=0;
0144 gloss=zeros(size(gg));
0145
0146 dginit=gg'*dd;
0147
0148 if dginit>=0
0149 if opt.display>0
0150 fprintf('dg=%g is not a descending direction!\n', dginit);
0151 end
0152 step = 0;
0153 ret = -1;
0154 return;
0155 end
0156
0157 Ip=find(dd>0);
0158 In=find(dd<0);
0159 step=min([step, 0.999*min((xx(In)-ll(In))./(-dd(In))), 0.999*min((uu(Ip)-xx(Ip))./dd(Ip))]);
0160
0161
0162 xx0 = xx;
0163 f0 = fval;
0164 gg0 = gg;
0165 cc = 0;
0166
0167 if opt.display>2
0168 fprintf('finit=%.20f\n',f0);
0169 end
0170
0171 while cc<opt.max_linesearch
0172 ftest = f0 + opt.ftol*step*dginit;
0173 xx = xx0 + step*dd;
0174
0175 if ~isempty(Ac)
0176 bineq = all(Ac*xx<=bc);
0177 else
0178 bineq = true;
0179 end
0180
0181 if bineq && all(xx>=ll) && all(xx<=uu)
0182 [fval, gg, info]=feval(fun, xx, info, varargin{:});
0183
0184 if fval<=ftest
0185 break;
0186 end
0187 else
0188 fval = inf;
0189 end
0190 if opt.display>2
0191 fprintf('[%d] step=%g fval=%.20f > ftest=%.20f\n', cc, step, fval, ftest);
0192 end
0193
0194 step = step/2;
0195 cc = cc+1;
0196 end
0197
0198 if cc==opt.max_linesearch
0199 if opt.display>0
0200 fprintf('Maximum linesearch=%d reached\n', cc);
0201 end
0202 xx = xx0;
0203 [fval, gg, info]=feval(fun, xx, info, varargin{:});
0204 step = 0;
0205 ret = -2;
0206 return;
0207 end
0208
0209
0210 ret = 0;