function [xk, xklist, alphaklist] = newtonbt(x0,f,tol)
% NEWTONBT  Newton optimization with back-tracking and an absolute
% tolerance on the norm of the gradient.  This algorithm is implied in
% section 2.2 of Nocedal & Wright.  If xk is the current iterate then the
% search vector is
%    pk = - H f (xk) \ grad f(xk).
% Back-tracking line search is applied to get alphak (bt.m = Algorithm 3.1).
% The next iterate is
%    xk <-- xk + alphak pk
%
% Usage:
%    [xk, xklist, alphaklist] = newtonbt(x0,f,tol)
% with inputs
%    x0          n-entry column vector with initial iterate (location)
%    f           function handle (e.g. "f" if f is anonymous or "@f" if f is
%                f.m); f() needs to return f(xk), gradient df(xk), and Hessian
%                Hf(xk) as first three outputs:    [fxk, dfxk, Hfxk] = f(xk)
%    tol         stop when norm of gradient is less than this number
% and outputs
%    xk          Nth iterate
%    xklist      all iterates as N+1 column matrix
%    alphaklist  step lengths as row vector of length N
%
% Example, as visualized in slides http://bueler.github.io/M661F16/pits.pdf:
%   x0 = [1.5 0.5]';
%   x = newtonbt(x0,@pits,1.0e-2)
%
% Example for Exercise 3.1 in Nocedal & Wright:
%   x0 = [1.2 1.2]';
%   [x, xlist, alphalist] = newtonbt(x0,@rosenbrock,1.0e-4)
%
% Requires: BT

maxiters = 100;   % never take more steps than this
xk = x0(:);       % force as column
xklist = [xk];
alphaklist = [];
for k = 1:maxiters
    [fxk, dfxk, Hfxk] = f(xk);
    if norm(dfxk) < tol          % absolute tolerance on gradient f
        break
    end
    pk = - Hfxk \ dfxk;          % nontrivial Newton step; \ is Gauss elim
    if dfxk' * pk < 0.0
         alphak = bt(xk,pk,f,dfxk);
    else
         warning('not descent direction ... unable to do backtracking')
         alphak = 1.0;
    end
    xk = xk + alphak * pk;
    xklist = [xklist xk];        % append to list
    alphaklist = [alphaklist alphak];
end
