function [x,fval,exitflag,output,grad]=lbfgs(funfcn,x_init,optim)
%FMINLBFGS finds a local minimum of a function of several variables.
% This optimizer is developed for image registration methods with large
% amounts of unknown variables.
%
% Optimization methods supported:
% - Quasi Newton BroydenFletcherGoldfarbShanno (BFGS)
% - Limited memory BFGS (L-BFGS)
% - Steepest Gradient Descent optimization.
%
% [X,FVAL,EXITFLAG,OUTPUT,GRAD] = FMINLBFGS(FUN,X0,OPTIONS)
%
% Inputs,
% FUN: Function handle or string which is minimized, returning an
% error value and optional the error gradient.
% X0: Initial values of unknowns can be a scalar, vector or matrix
% (optional)
% OPTIONS: Structure with optimizer options, made by a struct or
% optimset. (optimset doesnot support all input options)
%
% Outputs,
% X : The found location (values) which minimize the function.
% FVAL : The minimum found
% EXITFLAG : Gives value, which explain why the minimizer stopt
% OUTPUT : Structure with all important ouput values and parameters
% GRAD : The gradient at this location
%
% Extended description of input/ouput variables
% OPTIONS,
% OPTIONS.GoalsExactAchieve : If set to 0, a line search method is
% used which uses a few function calls to do a good line
% search. When set to 1 a normal line search method with Wolfe
% conditions is used (default).
% OPTIONS.GradConstr, Set this variable to true if gradient calls are
% cpu-expensive (default). If false more gradient calls are
% used and less function calls.
% OPTIONS.HessUpdate : If set to 'bfgs', BroydenFletcherGoldfarbShanno
% optimization is used (default), when the number of unknowns is
% larger then 3000 the function will switch to Limited memory BFGS,
% or if you set it to 'lbfgs'. When set to 'steepdesc', steepest
% decent optimization is used.
% OPTIONS.StoreN : Number of iterations used to approximate the Hessian,
% in L-BFGS, 20 is default. A lower value may work better with
% non smooth functions, because than the Hessian is only valid for
% a specific position. A higher value is recommend with quadratic equations.
% OPTIONS.GradObj : Set to 'on' if gradient available otherwise finited difference
% is used.
% OPTIONS.Display : Level of display. 'off' displays no output; 'plot' displays
% all linesearch results in figures. 'iter' displays output at each
% iteration; 'final' displays just the final output; 'notify'
% displays output only if the function does not converge;
% OPTIONS.TolX : Termination tolerance on x, default 1e-6.
% OPTIONS.TolFun : Termination tolerance on the function value, default 1e-6.
% OPTIONS.MaxIter : Maximum number of iterations allowed, default 400.
% OPTIONS.MaxFunEvals : Maximum number of function evaluations allowed,
% default 100 times the amount of unknowns.
% OPTIONS.DiffMaxChange : Maximum stepsize used for finite difference gradients.
% OPTIONS.DiffMinChange : Minimum stepsize used for finite difference gradients.
% OPTIONS.OutputFcn : User-defined function that an optimization function calls
% at each iteration.
% OPTIONS.rho : Wolfe condition on gradient (c1 on wikipedia), default 0.01.
% OPTIONS.sigma : Wolfe condition on gradient (c2 on wikipedia), default 0.9.
% OPTIONS.tau1 : Bracket expansion if stepsize becomes larger, default 3.
% OPTIONS.tau2 : Left bracket reduction used in section phase,
% default 0.1.
% OPTIONS.tau3 : Right bracket reduction used in section phase, default 0.5.
% FUN,
% The speed of this optimizer can be improved by also providing
% the gradient at X. Write the FUN function as follows
% function [f,g]=FUN(X)
% f , value calculation at X;
% if ( nargout > 1 )
% g , gradient calculation at X;
% end
% EXITFLAG,
% Possible values of EXITFLAG, and the corresponding exit conditions
% are
% 1, 'Change in the objective function value was less than the specified tolerance TolFun.';
% 2, 'Change in x was smaller than the specified tolerance TolX.';
% 3, 'Magnitude of gradient smaller than the specified tolerance';
% 4, 'Boundary fminimum reached.';
% 0, 'Number of iterations exceeded options.MaxIter or number of function evaluations exceeded options.FunEvals.';
% -1, 'Algorithm was terminated by the output function.';
% -2, 'Line search cannot find an acceptable point along the current search';
%
% Examples
% options = optimset('GradObj','on');
% X = fminlbfgs(@myfun,2,options)
%
% % where myfun is a MATLAB function such as:
% function [f,g] = myfun(x)
% f = sin(x) + 3;
% if ( nargout > 1 ), g = cos(x); end
%
% See also OPTIMSET, FMINSEARCH, FMINBND, FMINCON, FMINUNC, @, INLINE.
%
% Function is written by D.Kroon University of Twente (Updated Nov. 2010)
% Read Optimalisation Parameters
defaultopt = struct('Display','final','HessUpdate','bfgs','GoalsExactAchieve',1,'GradConstr',true, ...
'TolX',1e-6,'TolFun',1e-6,'GradObj','off','MaxIter',400,'MaxFunEvals',100*numel(x_init)-1, ...
'DiffMaxChange',1e-1,'DiffMinChange',1e-8,'OutputFcn',[], ...
'rho',0.0100,'sigma',0.900,'tau1',3,'tau2', 0.1, 'tau3', 0.5,'StoreN',20);
if (~exist('optim','var'))
optim=defaultopt;
else
f = fieldnames(defaultopt);
for i=1:length(f),
if (~isfield(optim,f{i})||(isempty(optim.(f{i})))), optim.(f{i})=defaultopt.(f{i}); end
end
end
% Initialize the data structure
data.fval=0;
data.gradient=0;
data.fOld=[];
data.xsizes=size(x_init);
data.numberOfVariables = numel(x_init);
data.xInitial = x_init(:);
data.alpha=1;
data.xOld=data.xInitial;
data.iteration=0;
data.funcCount=0;
data.gradCount=0;
data.exitflag=[];
data.nStored=0;
data.timeTotal=tic;
data.timeExtern=0;
% Switch to L-BFGS in case of more than 3000 unknown variables
if(optim.HessUpdate(1)=='b')
if(data.numberOfVariables<3000),
optim.HessUpdate='bfgs';
else
optim.HessUpdate='lbfgs';
end
end
if(optim.HessUpdate(1)=='l')
succes=false;
while(~succes)
try
data.deltaX=zeros(data.numberOfVariables,optim.StoreN);
data.deltaG=zeros(data.numberOfVariables,optim.StoreN);
data.saveD=zeros(data.numberOfVariables,optim.StoreN);
succes=true;
catch ME
warning('fminlbfgs:memory','Decreasing StoreN value because out of memory');
succes=false;
data.deltaX=[]; data.deltaG=[]; data.saveD=[];
optim.StoreN=optim.StoreN-1;
if(optim.StoreN<1)
rethrow(ME);
end
end
end
end
exitflag=[];
% Display column headers
if(strcmp(optim.Display,'iter'))
disp(' Iteration Func-count Grad-count f(x) Step-size');
end
% Calculate the initial error and gradient
data.initialStepLength=1;
[data,fval,grad]=gradient_function(data.xInitial,funfcn, data, optim);
data.gradient=grad;
data.dir = -data.gradient;
data.fInitial = fval;
data.fPrimeInitial= data.gradient'*data.dir(:);
data.fOld=data.fInitial;
data.xOld=data.xInitial;
data.gOld=data.gradient;
gNorm = norm(data.gradient,Inf); % Norm of gradient
data.initialStepLength = min(1/gNorm,5);
% Show the current iteration
if(strcmp(optim.Display,'iter'))
s=sprintf(' %5.0f %5.0f %5.0f %13.6g ',data.iteration,data.funcCount,data.gradCount,data.fInitial); disp(s);
end
% Hessian intialization
if(optim.HessUpdate(1)=='b')
data.Hessian=eye(data.numberOfVariables);
end
% Call output function
if(call_output_function(data,optim,'init')), exitflag=-1; end
% Start Minimizing
while(true)
% Update number of itterations
data.iteration=data.iteration+1;
