% [X,F] = newtoncg(f,paramf,x0[,tol,maxit,convcrit]) Optimisation by Newton-CG
%
% Reference: p. 136ff, 139ff in:
%   Nocedal and Wright: "Numerical Optimization", Springer, 1999.
%
% Example usage: see steepdesc.m.
%
% 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.
%   convcrit: convergence criterion for inexact Newton steps (p. 136ff);
%      one of 'linear', 'superlinear' or 'quadratic'. Default: 'quadratic'.
% Out:
%   X: ? x n list of iterates.
%   F: ? x n 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] = newtoncg(f,paramf,x0,tol,maxit,convcrit)

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

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

% $$$ fprintf('CG iterations:');
while (norm(g) > tol) & (k < maxit)
  
  % Convergence tolerance for inexact Newton step with linear CG
  switch convcrit
   case 'superlinear'
    eta = min(0.5,sqrt(norm(g)));
   case 'quadratic'
    eta = min(0.5,norm(g));
   otherwise
    eta = 0.5;			% 'linear'
  end
  
  % New search direction: inexact Newton step with linear CG
  [p,cgit] = lincg1(A,-g',zeros(size(x0)),eta*norm(g));
% $$$   fprintf(' %d',cgit);

% $$$   % Backtracking line search
% $$$   alpha0 = 1; rho = 0.8; c = 1e-4;	% Backtracking l. s. parameters
% $$$   alpha = alpha0;
% $$$   tmp = c*g*p;
% $$$   while f((x+alpha*p)',paramf{:}) > ff + alpha*tmp
% $$$     alpha = rho*alpha;
% $$$   end

  % No line search
  alpha = 1;
  
  x = x + alpha*p;		% New iterate
  [ff g A] = f(x',paramf{:});	% New value of function, gradient, Hessian
  X = [X;x']; F = [F;ff];	% Update lists
  k = k + 1;
  
end


% x = lincg1(A,b,x0[,tol,maxit])
% Solution of a linear system by conjugate gradients
%
% This is an adaptation of "lincg" to deal with non-positive-definite
% linear systems (as can arise in a Newton iteration).
%
% Reference: p. 111ff in:
%   Nocedal and Wright: "Numerical Optimization", Springer, 1999.
%
% In:
%   A: n x n positive definite matrix.
%   b: n x 1 vector.
%   x0: n x 1 vector, the starting point.
%   tol: tolerance to stop iterating (minimum gradient norm). Default: 1e-6.
%   maxit: maximum number of iterations. Default: n.
% Out:
%   x: n x 1 (approximate) solution.
%
% 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,k] = lincg1(A,b,x0,tol,maxit)

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

x = x0;				% Current iterate
r = A*x-b;			% Current gradient of the quadratic function
p = -r;				% Current search direction
k = 0;				% Number of iterations
tmp = 1;			% Curvature indicator

while (norm(r) > tol) & (k < maxit) & (tmp > 0)
  
  tmp = p'*A*p;
  if tmp > 0			% Positive curvature
    alpha = (r'*r)/tmp;
    x = x + alpha*p;		% New iterate
    oldr = r;
    r = r + alpha*A*p;		% New gradient
    beta = (r'*r)/(oldr'*oldr);
    p = -r + beta*p;		% New search direction
    k = k + 1;
  elseif k == 0			% Non-positive curvature at first iteration
    x = -r;
  end
  
end

