# HG changeset patch # User Jaroslav Hajek # Date 1237276148 -3600 # Node ID 22c8272af34bfcc66b1dc9be25b8cd6269812cd8 # Parent 193804a4f82f14a492daf2fe97c6bcffdac12475 improvements to fsolve & family diff -r 193804a4f82f -r 22c8272af34b scripts/ChangeLog --- a/scripts/ChangeLog Sun Mar 15 10:09:44 2009 +0100 +++ b/scripts/ChangeLog Tue Mar 17 08:49:08 2009 +0100 @@ -1,3 +1,13 @@ +2009-03-17 Jaroslav Hajek + + * 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 * statistics/base/var.m: a -> x. diff -r 193804a4f82f -r 22c8272af34b scripts/optimization/__dogleg__.m --- a/scripts/optimization/__dogleg__.m Sun Mar 15 10:09:44 2009 +0100 +++ b/scripts/optimization/__dogleg__.m Tue Mar 17 08:49:08 2009 +0100 @@ -53,6 +53,7 @@ alpha = 0; endif else + alpha = delta / xn; snm = 0; endif ## Form the appropriate convex combination. diff -r 193804a4f82f -r 22c8272af34b scripts/optimization/__fdjac__.m --- a/scripts/optimization/__fdjac__.m Sun Mar 15 10:09:44 2009 +0100 +++ b/scripts/optimization/__fdjac__.m Tue Mar 17 08:49:08 2009 +0100 @@ -17,21 +17,19 @@ ## . ## -*- 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 diff -r 193804a4f82f -r 22c8272af34b scripts/optimization/fsolve.m --- a/scripts/optimization/fsolve.m Sun Mar 15 10:09:44 2009 +0100 +++ b/scripts/optimization/fsolve.m Tue Mar 17 08:49:08 2009 +0100 @@ -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