gistfile1.matlab

function a = lsa(fun,x,d,a1,varargin)
% a = lsa(fun,x,d,a1)
%
% Implementation of Line Search Algorithm with Strong Wolfe conditions
% as found J. Nocedal, S. Wright, Numerical Optimization, 1999 edition
% Algorithm 3.2 on page 59
% 
% Output arguments:
%       a  : final stepsize
%
% Input arguments:
%       fun  : function handle as [f,g] = fun(x)
%       x    : point in which the line search is executed
%       d    : search direction
%       a1   : first step-length tried
%
% Optional input arguments (if one of these is provided but not the
% previous ones, use [] as a placeholder the previous):
%       c1   : sufficient decrease check (Wolfe1), default at '1e-4'
%       c2   : curvature condition (Wolfe2), default at '0.5'
%       rho  : to compute next trial step, default at '0.8'
%       amax : maximum trial step-length, default at '10*a1'
%       maxit: maximum number of trials, default at '100'
%
% author: Davide Taviani (23/06/2011)
% website: http://davidetaviani.com


% a1 = a_i
% a0 = a_{i-1}

% initialization of default parameters
c1 = 1e-4;
c2 = 0.5;
rho = 0.8;
amax = 10*a1;
maxit = 100;


if nargin > 3
    
    if ~isempty(varargin{1});
        c1 = varargin{1};
    end
    
    if nargin > 4
        % check on r
        if ~isempty(varargin{2});
            c2 = varargin{2};
        end
        
        if nargin > 5
            % check on c
            if ~isempty(varargin{3});
                rho = varargin{3};
                
            end
            
            if nargin > 6
                % check on a0
                if ~isempty(varargin{4});
                    amax = varargin{4};
                end
                
                if nargin > 7
                    if ~isempty(varargin{5})
                        maxit = varargin{5};
                    end
                end
            end
        end
    end
end


a0 = 0;                 % 0th steplength is 0
i=1;
[f0,g0] = fun(x);       % storing information of function and gradient
                        % at the point for further use

fold = fun(x+a0*d);     % initialization for function with the previous step-length

while 1
    [f,g] = fun(x+a1*d);    % function with current step-length
    if (f > f0+c1*a1*g0) | ((i>1) & f > fold)  % (sufficent decrease check or comparison between f and fold
        a = zoom_w(fun,x,d,a0,a1,c1,c2);  % a suitable step-length is in [a0,a1]
        return;
    end
    if abs(g) <= -c2*g0  % curvature condition
        a = a1;  % current step-length satisfies strong wolfe conditions
        return;
    end
    if g >= 0;   
        a = zoom_w(fun,x,d,a1,a0,c1,c2); % suitable step-length is in [a1,a0]
        return;
    end
    
    if i == maxit
        disp('Maximum number of iteration for Line Search reached');
        a = a1;
        return;
    end
    
    % update for next loop
    i=i+1;
    a0 = a1;
    a1 = rho*a0+(1-rho)*amax;
    fold = f;
        
end

gistfile2.matlab

function alpha = zoom_w(fun,x,d,alo,ahi,varargin)

% a = lsa(fun,x,d,a1)
%
% Implementation of Zoom algorithm 
% as found J. Nocedal, S. Wright, Numerical Optimization, 1999 edition
% Algorithm 3.3 on page 60
% 
% Output arguments:
%       alpha  : final stepsize
%
% Input arguments:
%       fun  : function handle as [f,g] = fun(x)
%       x    : point in which the line search is executed
%       d    : search direction
%       alo  : lowest bound of interval
%       ahi  : highest bound of interval
%
% Optional input arguments (if one of these is provided but not the
% previous ones, use [] as a placeholder the previous):
%       c1   : sufficient decrease check (Wolfe1), default at '1e-4'
%       c2   : curvature condition (Wolfe2), default at '0.5'
%       maxit: maximum number of trials, default at '100'
%
% author: Davide Taviani (23/06/2011)
% website: http://davidetaviani.com

c1 = 1e-4;
c2 = 0.5;
maxit = 20;

if nargin > 5
    
    if ~isempty(varargin{1});
        c1 = varargin{1};
    end
    
    if nargin > 6
        % check on r
        if ~isempty(varargin{2});
            c2 = varargin{2};
        end
        
        if nargin > 7
            % check on c
            if ~isempty(varargin{3});
                maxit = varargin{3};
                
            end
        end
    end
end

[f0,g0] = fun(x);
j = 0;
while 1
    a = (alo+ahi)/2;       % trial step-length is the middle point of [alo,ahi]
    [f,g] = fun(x+a*d);
    if (f > f0 + c1*a*g0 | f >= fun(x+alo*d)) % sufficient decrease or comparison with alo
        ahi = a;      % narrows the interval between [alo,ahi] 
    else
        if abs(g) <= -c2*g0  % curvature condition
            alpha = a;       % a satisfies the strong wolfe conditions
            return;
        end
        if g*(ahi-alo) >= 0
            ahi = alo;
        end
        alo = a;             % the interval is now [a,alo]
    end
    
    
    if j==maxit
        alpha = a;           % escape condition
        return;
    end
    
    j = j+1;
end

gistfile3.matlab

function [x,i] = sdd(fun,x0,maxit,varargin)
% [x,i] = sdd(fun,x0,maxit,tol,r,c,a0)
%
% Implementation of Steepest Descent with Backtracking
% 
% Output arguments:
%       x  : final point of method
%       i  : number of iteration to achieve precision 'tol'
%
% Input arguments:
%       fun    : function handle as [f,g] = fun(x)
%       x0     : starting point
%       maxit  : maximum number of iterations
%
% Optional input arguments (if one of these is provided but not the previous ones, use [] as a placeholder the previous):
%       tol    : stopping criterion on the norm of gradient at current x
%                default at '1e-3'.
%       r      : backtracking factor, default at '0.5'
%       c      : parameter for sufficient decrease check, default at '1e-4'
%       a0     : initial trial step, default at '1'.
%       
%
% author: Davide Taviani (23/06/2011)
% website: http://davidetaviani.com

% initialization of default parameters
tol = 1e-3;
r = 0.5;
c = 1e-4;
a0 = 1;

if nargin > 3
    
    % check on tol
    if ~isempty(varargin{1});
        tol = varargin{1};
    end
    if nargin > 4
        % check on r
        if ~isempty(varargin{2});
            r = varargin{2};
        end
        
        if nargin > 5
            % check on c
            if ~isempty(varargin{3});
                c = varargin{3};
                
            end
            
            if nargin > 6
                % check on a0
                if ~isempty(varargin{4});
                    a0 = varargin{4};
                    
                end
            end
        end
    end
end

x = x0(:);                      % first point
[~,g] = fun(x);                 % gradient at first point
for i=1:maxit
    d = -g;                     % computation of discent direction
    a = bkt(fun,c,r,a0,x,d);    % backtracking
    
    x = x+a*d;                  % computation of next point
    [~,g] = fun(x);             % gradient at next point
    
    % stopping criterion
    if norm(g) < tol
        fprintf('Convergence reached!');
        break;
    end
    
    % warning if tol is not reached in maxit steps
    if i == maxit
        disp('Maximum number of iteration reached');
    end
end

gistfile4.matlab

function [x,i] = sdd_w(fun,x0,maxit,varargin)
% [x,i] = sdd_w(fun,x0,maxit,tol)
%
% Implementation of Steepest Descent with Strong Wolfe Conditions Line Search
% 
% Output arguments:
%       x  : final point of method
%       i  : number of iteration to achieve precision 'tol'
%
% Input arguments:
%       fun    : function handle as [f,g] = fun(x)
%       x0     : starting point
%       maxit  : maximum number of iterations
%
% Optional input arguments:
%       tol    : stopping criterion on the norm of gradient at current x
%                default at '1e-3'.
%
% author: Davide Taviani (21/06/2011)
% website: http://davidetaviani.com

if ~isempty(varargin)
    tol = varargin{1};
else
    tol = 1e-3;
end

x = x0(:);
[f,g] = fun(x);
for i=1:maxit
    d = -g;
    a = lsa(fun,x,d,1);
    x = x+a*d;
    [f,g] = fun(x);
    if norm(g) < tol
        break;
        fprintf('Convergence reached! x=%g',x);
    end
    if i == maxit
        disp('Maximum number of iteration reached');
    end
end

gistfile5.matlab

function [x,k] = bfgs_w(fun,x0,H,tol)
% [x,k] = bfgs_w(fun,x0,H,tol)
%
% Implementation of BFGS with Strong Wolfe Conditions Line Search
% 
% Output arguments:
%       x  : final point of method
%       k  : number of iteration to achieve precision 'tol'
%
% Input arguments:
%       fun: function handle as [f,g] = fun(x)
%       x0 : starting point
%       H  : initial approximation of Hessian
%       tol: tolerance (stopping criterion on norm of gradient at x)
%
% author: Davide Taviani (23/06/2011)
% website: http://davidetaviani.com

n = size(H);
x = x0(:);
k = 0;
[~,g] = fun(x);
while norm(g) > tol
    d = -H*g;              % search direction computation
    a = lsa(fun,x,d,1);    % line search with stronge wolfe conditions
    x1 = x+a*d;            % compute next point x_{k+1}
    [~,g1] = fun(x1);      % gradient at x_{k+1}
    
    s = x1-x;
    y = g1-g;              % number required by Hessian update
    rho = 1/(y'*s);
    
   
    H = (eye(n)-rho*s*y')*H*(eye(n)-rho*y*s')+rho*(s*s');  % update of the Hessian estimate
    
   
    k = k+1;
    g = g1;                % update of the variables for next loop
    x = x1;
end
    

gistfile6.matlab

function [x,i] = newt(fun,x0,maxit,tol,varargin)
% newt(fun,x0,maxit,tol,varargin)
%
% Implementation of Modified Newton with backtracking
% 
% Output arguments:
%       x  : final point of method
%       k  : number of iteration to achieve precision 'tol'
%
% Input arguments:
%       fun   : function handle as [f,g] = fun(x)
%       x0    : starting point
%       maxit : maximum number of iterations
%       tol   : tolerance (stopping criterion on norm of gradient at x)
%
% Optional input arguments (if one of these is provided but not the previous ones, use [] as a placeholder the previous):
%       r      : backtracking factor, default at '0.5'
%       c      : parameter for sufficient decrease check, default at '1e-4'
%       a0     : initial trial step, default at '1'.
%
% author: Davide Taviani (23/06/2011)
% website: http://davidetaviani.com

r = 0.5;
c = 1e-4;
a0 = 1;

if nargin > 3
    
    % check on tol
    if ~isempty(varargin{1});
        tol = varargin{1};
    end
    if nargin > 4
        % check on r
        if ~isempty(varargin{2});
            r = varargin{2};
        end
        
        if nargin > 5
            % check on c
            if ~isempty(varargin{3});
                c = varargin{3};
                
            end
            
            if nargin > 6
                % check on a0
                if ~isempty(varargin{4});
                    a0 = varargin{4};
                    
                end
            end
        end
    end
end

x = x0(:);
for i=1:maxit
    [~,g,H] = fun(x);                   % computation of gradient and Hessian
    H = H+sqrt(eps)*(eye(size(H)));     % perturbation of the Hessian
    p = -H\g;                           % search direction computation
    a = bkt(fun,c,r,a0,x,p);            % backtracking
    x = x+a*p;                          % next point
    [~,g,H] = fun(x);                   % gradient-hessian at next point
    
    if norm(g) < tol                    % stopping criterion
        break;
    end
    
    if i==maxit                         % warning if tol is not reached in maxit steps
        disp('Maximum number of iteration reached');
    end
end