1;

##
## Demonstration of steepest ascent.
##

function y = fun(x, A)
  y = x'*A*x;
endfunction


function G = grad(x, A)
  G = 2*A*x;
endfunction


function x = opt_SA(x, A)

  ## The error from the previous iteration.
  f_lag = 0;

  while (1)

    G = grad(x, A);
    ng = norm(G);
    if (ng < 1e-6)
      break;
    endif

    ## The inner bracket point is at the current iterate.
    c1 = 0;
    f1 = fun(x, A);

    ## Print out the current gradient, function value, and the log ratio
    ## of consecutive errors.
    printf("|G|=%15.10f  f(x)=%15.12f  log(|E(i)/E(i-1)|)=%7.4f\n", ...
	   ng, f1, log(abs(f1)) - log(abs(f_lag)));
    f_lag = f1;

    ## Find an uphill point to be the middle bracket point.
    c2 = 1;
    while (1)
      x2 = x + c2*G;
      f2 = fun(x2, A);
      if (f2 > f1)
	break;
      endif
      c2 = c2/2;
    endwhile

    ## Find the outer bracket point.
    c3 = 2*c2;
    while (1)
      x3 = x + c3*G;
      f3 = fun(x3, A);
      if (f3 < f2)
	break;
      endif
      c3 = 2*c3;
    endwhile

    ## Bisection.
    while (c3-c1 > 1e-8)

      if (c3-c2 > c2-c1)

	cc = (c2+c3)/2;
	xx = x + cc*G;
	ff = fun(xx, A);
	
	if (ff > f2)
	  c1=c2;
	  f1=f2;
	  c2=cc;
	  f2=ff;
	else
	  c3=cc;
	  f3=ff;
	endif

      else

	cc = (c1+c2)/2;
	xx = x + cc*G;
	ff = fun(xx, A);

	if (ff > f2)
	  c3=c2;
	  f3=f2;
	  c2=cc;
	  f2=ff;
	else
	  c1=cc;
	  f1=ff;
	endif

      endif

    endwhile
    
    ## Move to the optimal point.
    x = x + c2*G;

  endwhile

endfunction


A = randn(5,5);
A = -A*A';

x = randn(5,1);

x = opt_SA(x, A);