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.
% This version is based on the modified C code L-BFGS-B-C, and so has
% a slightly different calling syntax than previous versions.
%
% The minimization problem that it 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: [] )
% Results saved in columns 3 and higher of info.err variable
%
% Stephen Becker, srbecker@alumni.caltech.edu
% Feb 14, 2012
% Updated Feb 21 2015, Stephen Becker, stephen.becker@colorado.edu
narginchk(3, 4)
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.
% just in case 'x' will be modified in-place
% (Feb 2015 version of code, it should not be modified,
% but just-in-case, may as well leave this )
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);
if ispc
nbd = int32(nbd);
else
nbd = int64(nbd);
end
% 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, 0 );
% 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, 0 ); % may crash if < 0
% 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 (and set to Inf to suppress)
printEvery = setOpts( 'printEvery', 1 );
errFcn = setOpts( 'errFcn', [] );
iprint = setOpts('verbose',-1);
% <0 for no output, 0 for some, 1 for more, 99 for more, 100 for more
% I recommend you set this -1 and use the Matlab print features
% (e.g., set printEvery )
fcn_wrapper([], [], [], 0); % initialize persistent variables
callF_wrapped = @(x,varargin) fcn_wrapper( callF, errFcn, maxIts, ...
printEvery, x, varargin{:} );
% callF_wrapped = @(x,varargin)callF(x); % also valid, but simpler
% Call the mex file
[f,x,taskInteger,outer_count, k] = lbfgsb_wrapper( m, x, l, u, nbd, ...
callF_wrapped, factr, pgtol, ...
iprint, maxIts, maxTotalIts);
info.iterations = outer_count;
info.totalIterations = k;
info.lbfgs_message1 = findTaskString( taskInteger );
errHist = fcn_wrapper([], [], [], printEvery);
info.err = errHist;
end % end of main function
function [f,g] = fcn_wrapper( callF, errFcn, maxIts, printEvery, x, varargin )
persistent k history
if isempty(k), k = 1; end
if nargin==4
% reset persistent variables and return information
if ~isempty(history) && ~isempty(k)
if printEvery > 0
printFcn(k,history);
end
f = history(1:k,:);
end
history = [];
k = [];
return;
end
if isempty( history )
width = 0;
if iscell( errFcn ), width = length(errFcn);
elseif ~isempty(errFcn), width = 1; end
width = width + 2; % include fcn and norm(grad) as well
history = zeros( maxIts, width );
end
% Find function value and gradient:
[f,g] = callF(x);
if nargin > 5
outerIter = varargin{1}+1;
history(outerIter,1) = f;
history(outerIter,2) = norm(g,Inf); % g is not projected
if isa( errFcn, 'function_handle' )
history(outerIter,3) = errFcn(x);
elseif iscell( errFcn )
for j = 1:length(errFcn)
history(outer_count,j+2) = errFcn{j}(x);
end
end
if outerIter > k
% Display info from *previous* input
% Since this may be called several times before outerIter
% is actually updated
% fprintf('At iterate %5d, f(x)= %.2e, ||grad||_infty = %.2e [MATLAB]\n',...
% k,history(k,1),history(k,2) );
if (printEvery > 0) && ~isinf(printEvery) && ~mod(k,printEvery)
printFcn(k,history);
end
k = outerIter;
end
end
end
function printFcn(k,history)
fprintf('Iter %5d, f(x) = %2e, ||grad||_infty = %.2e', ...
k, history(k,1), history(k,2) );
for col = 3:size(history,2)
fprintf(', %.2e', history(k,col) );
end
fprintf('\n');
end
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
function str = findTaskString( taskInteger )
% See the #define statements in lbfgsb.h
switch taskInteger
case 209
str = 'ERROR: N .LE. 0';
case 210
str = 'ERROR: M .LE. 0';
case 211
str = 'ERROR: FACTR .LT. 0';
case 3
str = 'ABNORMAL_TERMINATION_IN_LNSRCH.';
case 4
str = 'RESTART_FROM_LNSRCH.';
case 21
str = 'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL.';
case 22
str = 'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH.';
case 31
str = 'STOP: CPU EXCEEDING THE TIME LIMIT.';
case 32
str = 'STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIM.';
case 33
str = 'STOP: THE PROJECTED GRADIENT IS SUFFICIENTLY SMALL.';
case 101
str = 'WARNING: ROUNDING ERRORS PREVENT PROGRESS';
case 102
str = 'WARNING: XTOL TEST SATISIED';
case 103
str = 'WARNING: STP = STPMAX';
case 104
str = 'WARNING: STP = STPMIN';
case 201
str = 'ERROR: STP .LT. STPMIN';
case 202
str = 'ERROR: STP .GT. STPMAX';
case 203
str = 'ERROR: INITIAL G .GE. ZERO ';
case 204
str = 'ERROR: FTOL .LT. ZERO';
case 205
str = 'ERROR: GTOL .LT. ZERO';
case 206
str = 'ERROR: XTOL .LT. ZERO';
case 207
str = 'ERROR: STPMIN .LT. ZERO';
case 208
str = 'ERROR: STPMAX .LT. STPMIN';
case 212
str = 'ERROR: INVALID NBD';
case 213
str = 'ERROR: NO FEASIBLE SOLUTION';
otherwise
str = 'UNRECOGNIZED EXIT FLAG';
end
end