image thumbnail

LBFGSB (L-BFGS-B) mex wrapper

by

 

16 Feb 2012 (Updated )

Mex wrapper for lbfgsb v3.0 fortan library. L-bfgs-b solves box-constrained optimization.

lbfgsb
function [x,f,info] = lbfgsb( fcn, l, u, opts )
% x = lbfgsb( fcn, l, u )
%   uses the lbfgsb v.3.0 library (fortran files must be installed;
%       see compile_mex.m ) which is the L-BFGS-B algorithm.
%   The algorithm is similar to the L-BFGS quasi-Newton algorithm,
%   but also handles bound constraints via an active-set type iteration.
%
%  The minimization problem that is solves is:
%       min_x  f(x)     subject to   l <= x <= u
%
% 'fcn' is a function handle that accepts an input, 'x',
%   and returns two outputs, 'f' (function value), and 'g' (function gradient).
%
% 'l' and 'u' are column-vectors of constraints. Set their values to Inf
%   if you want to ignore them. (You can set some values to Inf, but keep
%   others enforced).
%
% The full format of the function is:
% [x,f,info] = lbfgsb( fcn, l, u, opts )
%   where the output 'f' has the value of the function f at the final iterate
%   and 'info' is a structure with useful information
%       (self-explanatory, except for info.err. The first column of info.err
%        is the history of the function values f, and the second column
%        is the history of norm( gradient, Inf ).  )
%
%   The 'opts' structure allows you to pass further options.
%   Possible field name values:
%
%       opts.x0     The starting value (default: all zeros)
%       opts.m      Number of limited-memory vectors to use in the algorithm
%                       Try 3 <= m <= 20. (default: 5 )
%       opts.factr  Tolerance setting (see this source code for more info)
%                       (default: 1e7 ). This is later multiplied by machine epsilon
%       opts.pgtol  Another tolerance setting, relating to norm(gradient,Inf)
%                       (default: 1e-5)
%       opts.maxits         How many iterations to allow (default: 100)
%       opts.maxTotalIts    How many iterations to allow, including linesearch iterations
%                       (default: 5000)
%       opts.printEvery     How often to display information (default: 1)
%       opts.errFcn         A function handle (or cell array of several function handles)
%                       that computes whatever you want. The output will be printed
%                       to the screen every 'printEvery' iterations. (default: [] )
%       opts.outputFcn      Similar to 'errFcn', but will save the output to
%                       the info.err variable.
%
% Stephen Becker, srbecker@alumni.caltech.edu
% Feb 14, 2012




error(nargchk(3, 4, nargin, 'struct'))
if nargin < 4, opts = struct([]); end

% Matlab doesn't let you use the .name convention with structures
%   if they are empty, so in that case, make the structure non-empty:
if isempty(opts), opts=struct('a',1) ; end

function out = setOpts( field, default, mn, mx )
    if ~isfield( opts, field )
        opts.(field)    = default;
    end
    out = opts.(field);
    if nargin >= 3 && ~isempty(mn) && any(out < mn), error('Value is too small'); end
    if nargin >= 4 && ~isempty(mx) && any(out > mx), error('Value is too large'); end
    opts    = rmfield( opts, field ); % so we can do a check later
end

% [f,g] = callF( x );
if iscell(fcn)
    % the user has given us separate functions to compute
    %   f (function) and g (gradient)
    callF   = @(x) fminunc_wrapper(x,fcn{1},fcn{2} );
else
    callF   = fcn;
end


n   = length(l);
if length(u) ~= length(l), error('l and u must be same length'); end
x0  = setOpts( 'x0', zeros(n,1) );
x   = x0 + 0; % important: we want Matlab to make a copy of this.
              %  'x' will be modified in-place

if size(x0,2) ~= 1, error('x0 must be a column vector'); end
if size(l,2) ~= 1, error('l must be a column vector'); end
if size(u,2) ~= 1, error('u must be a column vector'); end
if size(x,1) ~= n, error('x0 and l have mismatchig sizes'); end
if size(u,1) ~= n, error('u and l have mismatchig sizes'); end

% Number of L-BFGS memory vectors
% From the fortran driver file:
% "Values of m < 3  are not recommended, and
%  large values of m can result in excessive computing time.
%  The range  3 <= m <= 20 is recommended.  "
m   = setOpts( 'm', 5, 0 );


% 'nbd' is 0 if no bounds, 1 if lower bound only,
%       2 if both upper and lower bounds, and 3 if upper bound only.
% This .m file assumes l=-Inf and u=+Inf imply that there are no constraints.
% So, convert this to the fortran convention:
nbd     = isfinite(l) + isfinite(u) + 2*isinf(l).*isfinite(u);


% Some scalar settings, "factr" and "pgtol"
% Their descriptions, from the fortran file:

%     factr is a DOUBLE PRECISION variable that must be set by the user.
%       It is a tolerance in the termination test for the algorithm.
%       The iteration will stop when
%
%        (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch
%
%       where epsmch is the machine precision which is automatically
%       generated by the code. Typical values for factr on a computer
%       with 15 digits of accuracy in double precision are:
%       factr=1.d+12 for low accuracy;
%             1.d+7  for moderate accuracy;
%             1.d+1  for extremely high accuracy.
%       The user can suppress this termination test by setting factr=0.
factr   = setOpts( 'factr', 1e7 );

%     pgtol is a double precision variable.
%       On entry pgtol >= 0 is specified by the user.  The iteration
%         will stop when
%
%                 max{|proj g_i | i = 1, ..., n} <= pgtol
%
%         where pg_i is the ith component of the projected gradient.
%       The user can suppress this termination test by setting pgtol=0.
pgtol   = setOpts( 'pgtol', 1e-5 );


% Maximum number of outer iterations
maxIts  = setOpts( 'maxIts', 100, 1 );

% Maximum number of total iterations
%   (this includes the line search steps )
maxTotalIts     = setOpts( 'maxTotalIts', 5e3 );

% Print out information this often:
printEvery  = setOpts( 'printEvery', 1 );

errFcn      = setOpts( 'errFcn', [] );
outputFcn   = setOpts( 'outputFcn', [] );
width       = 0;
if iscell( outputFcn ), width = length(outputFcn);
elseif ~isempty(outputFcn), width = 1; end
width       = width + 2; % include fcn and norm(grad) as well
err         = zeros( maxIts, width );

% Make the work arrays
wa      = ones(2*m*n + 5*n + 11*m*m + 8*m,1);
iwa     = ones(3*n,1,'int32');
task    = 'START';
iprint  = 0;
csave   = '';
lsave   = zeros(4,1);
isave   = zeros(44,1, 'int32');
dsave   = zeros(29,1);
f       = 0;
g       = zeros(n,1);


outer_count     = 0;
for k = 1:maxTotalIts

    % Call the mex file. The way it works is that you call it,
    %   then it returns a "task". If that task starts with 'FG',
    %   it means it is requesting you to compute the function and gradient,
    %   and then call the function again.
    % If it is 'NEW_X', it means it has completed one full iteration.
    [f, task, csave, lsave, isave, dsave] = ...
        lbfgsb_wrapper( m, x, l, u, nbd, f, g, factr, pgtol, wa, iwa, task,iprint,...
        csave, lsave, isave, dsave );

    task    = deblank(task(1:60)); % this is critical!
                                   %otherwise, fortran interprets the string incorrectly

    if 1 == strfind( task, 'FG' )
        % L-BFGS-B requests that we compute the gradient and function value

        [f,g] = callF( x );

    elseif 1 == strfind( task, 'NEW_X' )
        outer_count     = outer_count + 1;

        % Display information if requested
        if ~mod( outer_count, printEvery )
            fprintf('Iteration %4d, f = %5.2e, ||g||_inf = %5.2e', ...
                outer_count, f, norm(g,Inf) );
            if isa( errFcn, 'function_handle' )
                fprintf('; error %.2e', errFcn(x) );
            elseif iscell( errFcn )
                for j = 1:length(errFcn)
                    fprintf('; err %.2e', errFcn{j}(x) );
                end
            end
            fprintf('\n');
        end

        err(outer_count,1)  = f;
        err(outer_count,2)  = norm(g,Inf);


        % Record information for the output, if requested
        % e.g. outputFcn = errFcn
        if isa( outputFcn, 'function_handle' )
            err(outer_count,3) = outputFcn(x);
        elseif iscell( outputFcn )
            for j = 1:length(outputFcn)
                err(outer_count,k+2) = outputFcn{j}(x);
            end
        end


    else
        break;
    end
end
info.err    = err(1:outer_count,:);
info.iterations     = outer_count;
info.totalIterations = k;
info.lbfgs_message1  = task;
info.lbfgs_message2  = csave;
info.g  = g;

end % end of main function


function [f,g] = fminunc_wrapper(x,F,G)
% [f,g] = fminunc_wrapper( x, F, G )
%   for use with Matlab's "fminunc"
f = F(x);
if nargin > 2 && nargout > 1
    g = G(x);
end

end
Error using ==> lbfgsb at 52
Not enough input arguments.

Contact us