Home > dal > lbfgs.m

lbfgs

PURPOSE ^

lbfgs - L-BFGS algorithm

SYNOPSIS ^

function [xx, status] = lbfgs(fun, xx, ll, uu, Ac, bc, info, opt, varargin)

DESCRIPTION ^

 lbfgs - L-BFGS algorithm

 Syntax:
  [xx, status] = lbfgs(fun, xx, ll, uu, <opt>)

 Input:
  fun     - objective function
  xx      - Initial point for optimization
  ll      - lower bound on xx
  uu      - upper bound on xx
  Ac      - inequality constraint:
  bc      -     Ac*xx<=bc
  opt     - Struct or property/value list of optional properties:
   .m          - size of limited memory
   .epsg       - gradient tolerance
   .maxiter    - maximum number of iterations
   .display    - display level
 
 Output:
  xx      - Final point of optimization
  status  - Various numbers

 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 % lbfgs - L-BFGS algorithm
0002 %
0003 % Syntax:
0004 %  [xx, status] = lbfgs(fun, xx, ll, uu, <opt>)
0005 %
0006 % Input:
0007 %  fun     - objective function
0008 %  xx      - Initial point for optimization
0009 %  ll      - lower bound on xx
0010 %  uu      - upper bound on xx
0011 %  Ac      - inequality constraint:
0012 %  bc      -     Ac*xx<=bc
0013 %  opt     - Struct or property/value list of optional properties:
0014 %   .m          - size of limited memory
0015 %   .epsg       - gradient tolerance
0016 %   .maxiter    - maximum number of iterations
0017 %   .display    - display level
0018 %
0019 % Output:
0020 %  xx      - Final point of optimization
0021 %  status  - Various numbers
0022 %
0023 % Copyright(c) 2009 Ryota Tomioka
0024 % This software is distributed under the MIT license. See license.txt
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 % Limited memory
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 % The initial step is gradient
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   % Perform line search
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   % Progress report
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 % || gnorm<opt.epsg
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   % L-BFGS update
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   % Initially set the negative gradient as descent direction
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;

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