## Demonstrate the following optimization methods:
##
## 1. Steepest ascent. 
## 2. Cyclic line searches over the canonical basis.
## 3. Fletcher-Powell conjugate gradient method. 
## 4. Polak-Ribiere conjugate gradient method.
##
## using the objective function
##
## f(x) = -x'Ax + b'x + c/(1 + x'x)
##
## The final result is the array GN containing the log 10 gradient norms
## for each method, and the array FV containing the function values for
## each method.


## The dimension of x.
d = 10;

## The scale of the non-quadratic term.
c = 1;


## Generate a SPD coefficient matrix (so that the maximum exists).
A = randn(2*d,d);
A = A'*A;

## Generate the linear term.
b = rand(d,1);


## Calculate the objective function and gradient.
function [V,G] = fun(x, A, b, c)

  ssx = sumsq(x);
  V = -x'*A*x + b'*x + c/(1 + ssx);
  G = -2*A*x + b - 2*c*x/(1 + ssx)^2;

endfunction


## A starting value.
xi = rand(d,1);

## Do the optimization using four methods.
for M=1:4

  ## Each algorithm starts in the same place.
  x = xi;

  ## Initialize the gradient and search direction.
  [u,g] = fun(xi, A, b, c);
  h = g;

  ## Initialize the search direction for the canonical basis search. 
  if (M == 2)
    g = [1; zeros(d-1,1)];
    h = g;
  endif

  ## Iterations loop for each method.
  for it=1:2*d
    
    ##
    ## Line search in the direction of h.
    ##
    
    ## Get an initial bracket.
    L = [-1,0,1];
    F = [fun(x+L(1)*h, A, b, c), fun(x, A, b, c), fun(x+L(3)*h, A, b, c)];
    while (F(2) < max(F(1), F(3)))
      L = L*2;
      F(1) = fun(x + L(1)*h, A, b, c);
      F(3) = fun(x + L(3)*h, A, b, c);
    endwhile
    
    ## Get the angle between successive search directions.
    if (it > 1)
      SA(it,M) = dot(h, h_old) / (norm(h)*norm(h_old));
    endif

    ## Check for conjugacy.
    if (it > 1)
      QA(it,M) = h'*A*h_old / (norm(h)*norm(h_old));
    endif

    ## Bisection.
    while (L(3) - L(1) > 1e-8)
      
      if (L(3) - L(2) > L(2) - L(1))
 	L1 = (L(2) + L(3)) / 2;
 	F1 = fun(x + L1*h, A, b, c);
 	if (F1 > F(2))
 	  L = [L(2), L1, L(3)];
 	  F = [F(2), F1, F(3)];
 	else
 	  L = [L(1), L(2), L1];
 	  F = [F(1), F(2), F1];
 	endif
      else
 	L1 = (L(1) + L(2)) / 2;
 	F1 = fun(x + L1*h, A, b, c);
 	if (F1 > F(2))
 	  L = [L(1), L1, L(2)];
 	  F = [F(1), F1, F(2)];
 	else
 	  L = [L1, L(2), L(3)];
 	  F = [F1, F(2), F(3)];
 	endif
      endif
      
    endwhile

    ## Move the point.
    x = x + L(2)*h;
    
    ##
    ## Update the search direction.
    ##
    
    ## Save the current direction.
    h_old = h;

    ## Get the new gradient.
    [u, g1] = fun(x, A, b, c);
    
    ## Save the final function value and gradient norm.
    GN(it,M) = sumsq(g1);
    FV(it,M) = u;

    ## Stepest ascent.
    if (M == 1)
      h = g1;
      g = g1;
      
    ## Search along the coordinate basis.
    elseif (M == 2)
      ii = find(h == 1);
      j = rem(ii(1), d);
      g = zeros(d,1);
      g(j+1) = 1;
      h = g;
      
    ## Fletcher-Powell conjugate gradient.
    elseif (M == 3)
      gamma = sumsq(g1) / sumsq(g);
      h = g1 + gamma*h;
      g = g1;
      
    ## Polak-Ribiere conjugate gradient.
    elseif (M == 4)
      gamma = dot(g1-g, g1) / sumsq(g);
      h = g1 + gamma*h;
      g = g1;
      
    endif

  endfor
  
endfor

## Put the gradient norms on the log 10 scale.
GN = log(GN) / log(10);