% [X,F] = steepdesc(f,paramf,x0[,tol,maxit]) Optimisation by steepest descent
%
% Reference: p. 21ff, 47ff in:
%   Nocedal and Wright: "Numerical Optimization", Springer, 1999.
%
% Example usage: the following will optimise Rosenbrock's function (with
% default parameter a = 100) starting at point [-1.2;1], then plot the
% function contours and the iterates, then estimate the convergence rate:
%
%   [X,F] = steepdesc(@rosenbrock,{},[-1.2;1]);
%   figure(1); fcontours(@rosenbrock,{},{},{},{},{},[-2 2;-1 3]); plotseq(X);
%   figure(2); [order,M] = convseq(X,[1 1]);
%
% Ditto for a quadratic function:
%
%   A = [1 0.1;0.1 2]; b = [-3;2]; c = 1;
%   [X,F] = steepdesc(@quadf,{A,b,c},[7;6]);
%   figure(1); fcontours(@quadf,{A,b,c},{},{},{},{},[-8 8;-8 8]); plotseq(X);
%   figure(2); fcontours(@quadf,{A,b,c},{},{},{},{},[-8 8;-8 8],F); plotseq(X);
%   figure(3); [order,M] = convseq(X,-(A \ b)'); cond(A)
%
% In:
%   f: a handle to a function that takes as input a list of vectors and
%      returns 3 lists: the function values, the gradients and the Hessians
%      (see rosenbrock.m for an example in 2D).
%   paramf: cell array containing the parameters for f. Use {} to pass no
%      parameters, or to use f's default parameters.
%   x0: n x 1 vector, the starting point.
%   tol: tolerance to stop iterating (minimum gradient norm). Default: 1e-5.
%   maxit: maximum number of iterations. Default: 1000.
% Out:
%   X: ? x n list of iterates.
%   F: ? x 1 list of the function value at the iterates.
%
% Any non-mandatory argument can be given the value [] to force it to take
% its default value.

% Copyright (c) 2005 by Miguel A. Carreira-Perpinan

function [X,F] = steepdesc(f,paramf,x0,tol,maxit)

% ---------- Argument defaults ----------
if ~exist('tol','var') | isempty(tol) tol = 1e-5; end;
if ~exist('maxit','var') | isempty(maxit) maxit = 1000; end;
% ---------- End of "argument defaults" ----------

% Backtracking line search parameters
alpha0 = 1; rho = 0.8; c = 1e-4;

x = x0;				% Current iterate
[ff g] = f(x',paramf{:});	% Current value of function and gradient
p = -g';			% Current search direction
X = x'; F = ff;			% Lists of iterates and their function values
k = 0;				% Number of iterations

while (norm(g) > tol) & (k < maxit)
  
  % Backtracking line search (see procedure 3.1 in the book):
  alpha = alpha0;
  tmp = c*g*p;
  while f((x+alpha*p)',paramf{:}) > ff + alpha*tmp
    alpha = rho*alpha;
  end
  
  x = x + alpha*p;		% New iterate
  [ff g] = f(x',paramf{:});	% New value of function and gradient
  p = -g';			% New search direction
  X = [X;x']; F = [F;ff];	% Update lists
  k = k + 1;
  
end

