--- a/scripts/ChangeLog
+++ b/scripts/ChangeLog
@@ -1,3 +1,13 @@
+2009-03-17  Jaroslav Hajek  <highegg@gmail.com>
+
+	* optimization/__fdjac__.m: Pass in fvec to save one evaluation.
+	* optimization/fsolve.m: Avoid redundant reevaluation when using
+	FD jacobians. Document how it can be done with user jacobians.  Make
+	first iteration special and call outputfcn after it. Skip updates
+	unless two successful iterations have occured.
+	* optimization/__dogleg__.m: Add missing alpha in the zero-gradient
+	case.
+
 2009-03-14  Jaroslav Hajek  <highegg@gmail.com>
 
 	* statistics/base/var.m: a -> x.

--- a/scripts/optimization/__dogleg__.m
+++ b/scripts/optimization/__dogleg__.m
@@ -53,6 +53,7 @@
         alpha = 0;
       endif
     else
+      alpha = delta / xn;
       snm = 0;
     endif
     ## Form the appropriate convex combination.

--- a/scripts/optimization/__fdjac__.m
+++ b/scripts/optimization/__fdjac__.m
@@ -17,21 +17,19 @@
 ## <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
-## @deftypefn{Function File} {[@var{fval}, @var{fjac}]} =  __fdjac__ (@var{fcn}, @var{x}, @var{err})
+## @deftypefn{Function File} {} __fdjac__ (@var{fcn}, @var{x}, @var{fvec}, @var{err})
 ## Undocumented internal function.
 ## @end deftypefn
 
-function [fvec, fjac] = __fdjac__ (fcn, x, err = 0)
+function fjac = __fdjac__ (fcn, x, fvec, err = 0)
   err = sqrt (max (eps, err));
-  fvec = fcn (x);
-  fv = fvec(:);
   h = abs (x(:))*err;
   h(h == 0) = err;
-  fjac = zeros (length (fv), numel (x));
+  fjac = zeros (length (fvec), numel (x));
   for i = 1:numel (x)
     x1 = x;
     x1(i) += h(i);
-    fjac(:,i) = (fcn (x1)(:) - fv) / h(i);
+    fjac(:,i) = (fcn (x1)(:) - fvec) / h(i);
   endfor
 endfunction
 

--- a/scripts/optimization/fsolve.m
+++ b/scripts/optimization/fsolve.m
@@ -72,10 +72,48 @@
 ## @item -3
 ## The trust region radius became excessively small. 
 ## @end table
-##
+## 
 ## Note: If you only have a single nonlinear equation of one variable, using 
 ## @code{fzero} is usually a much better idea.
 ## @seealso{fzero, optimset}
+##
+## Note about user-supplied jacobians:
+## As an inherent property of the algorithm, jacobian is always requested for a
+## solution vector whose residual vector is already known, and it is the last
+## accepted successful step. Often this will be one of the last two calls, but
+## not always. If the savings by reusing intermediate results from residual
+## calculation in jacobian calculation are significant, the best strategy is to
+## employ OutputFcn: After a vector is evaluated for residuals, if OutputFcn is
+## called with that vector, then the intermediate results should be saved for
+## future jacobian evaluation, and should be kept until a jacobian evaluation
+## is requested or until outputfcn is called with a different vector, in which
+## case they should be dropped in favor of this most recent vector. A short
+## example how this can be achieved follows:
+## @example
+## function [fvec, fjac] = my_optim_func (x, optimvalues, state)
+## persistent sav = [], sav0 = [];
+## if (nargin == 1)
+##   ## evaluation call
+##   if (nargout == 1)
+##     sav0.x = x; # mark saved vector
+##     ## calculate fvec, save results to sav0.
+##   elseif (nargout == 2)
+##     ## calculate fjac using sav.
+##   endif
+## else
+##   ## outputfcn call.
+##   if (all (x == sav0.x))
+##     sav = sav0;
+##   endif
+##   ## maybe output iteration status etc.
+## endif
+## endfunction
+##
+## ## ....
+## 
+## fsolve (@my_optim_func, x0, optimset ("OutputFcn", @my_optim_func, @dots{}))
+## @end example
+###
 ## @end deftypefn
 
 ## PKG_ADD: __all_opts__ ("fsolve");
@@ -129,11 +167,34 @@
   autodg = true;
 
   niter = 1;
-  nfev = 0;
+  nfev = 1;
 
   x = x0(:);
   info = 0;
 
+  ## Initial evaluation.
+  fvec = fcn (reshape (x, xsiz));
+  fsiz = size (fvec);
+  fvec = fvec(:);
+  fn = norm (fvec);
+  m = length (fvec);
+  n = length (x);
+
+  if (! isempty (outfcn))
+    optimvalues.iter = niter;
+    optimvalues.funccount = nfev;
+    optimvalues.fval = fn;
+    optimvalues.searchdirection = zeros (n, 1);
+    state = 'init';
+    stop = outfcn (x, optimvalues, state);
+    if (stop)
+      info = -1;
+      break;
+    endif
+  endif
+
+  nsuciter = 0;
+
   ## Outer loop.
   while (niter < maxiter && nfev < maxfev && ! info)
 
@@ -145,16 +206,12 @@
       if (issparse (fjac))
         updating = false;
       endif
+      fvec = fvec(:);
       nfev ++;
     else
-      [fvec, fjac] = __fdjac__ (fcn, reshape (x, xsiz));
-      nfev += 1 + length (x);
+      fjac = __fdjac__ (fcn, reshape (x, xsiz), fvec);
+      nfev += length (x);
     endif
-    fsiz = size (fvec);
-    fvec = fvec(:);
-    fn = norm (fvec);
-    m = length (fvec);
-    n = length (x);
 
     ## For square and overdetermined systems, we update a QR
     ## factorization of the jacobian to avoid solving a full system in each
@@ -199,6 +256,7 @@
       ## of course, how to define the above terms :)
     endif
 
+    lastratio = 0;
     nfail = 0;
     nsuc = 0;
     ## Inner loop.
@@ -249,7 +307,7 @@
       endif
 
       ## Update delta.
-      if (ratio < 0.1)
+      if (ratio < min(max(0.1, lastratio), 0.9))
         nsuc = 0;
         nfail ++;
         delta *= 0.5;
@@ -259,6 +317,7 @@
           break;
         endif
       else
+        lastratio = ratio;
         nfail = 0;
         nsuc ++;
         if (abs (1-ratio) <= 0.1)
@@ -274,6 +333,7 @@
         xn = norm (dg .* x);
         fvec = fvec1;
         fn = fn1;
+        nsuciter ++;
       endif
 
       niter ++;
@@ -321,7 +381,7 @@
       endif
 
       ## Criterion for recalculating jacobian.
-      if (! updating || nfail == 2)
+      if (! updating || nfail == 2 || nsuciter < 2)
         break;
       endif
 
@@ -347,6 +407,7 @@
   fvec = reshape (fvec, fsiz);
 
   output.iterations = niter;
+  output.successful = nsuciter;
   output.funcCount = nfev;
 
 endfunction