Mercurial > hg > octave-nkf
changeset 9212:6feb27c38da1
support central differences in fminunc and fsolve
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Mon, 18 May 2009 09:46:49 +0200 |
parents | f0c3d3fc4903 |
children | b9986ce0047c |
files | scripts/ChangeLog scripts/optimization/__fdjac__.m scripts/optimization/fminunc.m scripts/optimization/fsolve.m |
diffstat | 4 files changed, 38 insertions(+), 18 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,10 @@ +2009-05-14 Jaroslav Hajek <highegg@gmail.com> + + * optimization/__fdjac__.m: Support central differences. + * optimization/fsolve.m: Support central differences. Add FinDiffType + option. + * optimization/fminunc.m: Ditto. + 2009-05-17 Rik Wehbring <rdrider0-list@yahoo.com> * *.m: Simplify Texinfo documentation in .m scripts by removing
--- a/scripts/optimization/__fdjac__.m +++ b/scripts/optimization/__fdjac__.m @@ -21,16 +21,27 @@ ## Undocumented internal function. ## @end deftypefn -function fjac = __fdjac__ (fcn, x, fvec, err = 0) - err = sqrt (max (eps, err)); - h = abs (x(:))*err; - h(h == 0) = err; - fjac = zeros (length (fvec), numel (x)); - for i = 1:numel (x) - x1 = x; - x1(i) += h(i); - fjac(:,i) = (fcn (x1)(:) - fvec) / h(i); - endfor +function fjac = __fdjac__ (fcn, x, fvec, cdif, err = 0) + if (cdif) + err = (max (eps, err)) ^ (1/3); + h = max (abs (x), 1)*err; # FIXME? + fjac = zeros (length (fvec), numel (x)); + for i = 1:numel (x) + x1 = x2 = x; + x1(i) += h(i); + x2(i) -= h(i); + fjac(:,i) = (fcn (x1)(:) - fcn (x2)(:)) / (x1(i) - x2(i)); + endfor + else + err = sqrt (max (eps, err)); + h = max (abs (x), 1)*err; # FIXME? + fjac = zeros (length (fvec), numel (x)); + for i = 1:numel (x) + x1 = x; + x1(i) += h(i); + fjac(:,i) = (fcn (x1)(:) - fvec) / (x1(i) - x(i)); + endfor + endif endfunction
--- a/scripts/optimization/fminunc.m +++ b/scripts/optimization/fminunc.m @@ -29,10 +29,10 @@ ## @var{x0} determines a starting guess. The shape of @var{x0} is preserved ## in all calls to @var{fcn}, but otherwise it is treated as a column vector. ## @var{options} is a structure specifying additional options. -## Currently, @code{fsolve} recognizes these options: +## Currently, @code{fminunc} recognizes these options: ## @code{"FunValCheck"}, @code{"OutputFcn"}, @code{"TolX"}, ## @code{"TolFun"}, @code{"MaxIter"}, @code{"MaxFunEvals"}, -## @code{"GradObj"}. +## @code{"GradObj"}, @code{"FinDiffType"}. ## ## If @code{"GradObj"} is @code{"on"}, it specifies that @var{fcn}, ## called with 2 output arguments, also returns the Jacobian matrix @@ -74,7 +74,7 @@ x = optimset ("MaxIter", 400, "MaxFunEvals", Inf, \ "GradObj", "off", "TolX", 1.5e-8, "TolFun", 1.5e-8, "OutputFcn", [], "FunValCheck", "off", - "ComplexEqn", "off"); + "FinDiffType", "central"); return; endif @@ -90,6 +90,7 @@ n = numel (x0); has_grad = strcmpi (optimget (options, "GradObj", "off"), "on"); + cdif = strcmpi (optimget (options, "FinDiffType", "central"), "central"); maxiter = optimget (options, "MaxIter", 400); maxfev = optimget (options, "MaxFunEvals", Inf); outfcn = optimget (options, "OutputFcn"); @@ -152,8 +153,8 @@ grad = grad(:); nfev ++; else - grad = __fdjac__ (fcn, reshape (x, xsiz), fval)(:); - nfev += length (x); + grad = __fdjac__ (fcn, reshape (x, xsiz), fval, cdif)(:); + nfev += (1 + cdif) * length (x); endif if (niter == 1)
--- a/scripts/optimization/fsolve.m +++ b/scripts/optimization/fsolve.m @@ -126,7 +126,7 @@ x = optimset ("MaxIter", 400, "MaxFunEvals", Inf, \ "Jacobian", "off", "TolX", 1.5e-8, "TolFun", 1.5e-8, "OutputFcn", [], "Updating", "on", "FunValCheck", "off", - "ComplexEqn", "off"); + "ComplexEqn", "off", "FinDiffType", "central"); return; endif @@ -142,6 +142,7 @@ n = numel (x0); has_jac = strcmpi (optimget (options, "Jacobian", "off"), "on"); + cdif = strcmpi (optimget (options, "FinDiffType", "central"), "central"); maxiter = optimget (options, "MaxIter", 400); maxfev = optimget (options, "MaxFunEvals", Inf); outfcn = optimget (options, "OutputFcn"); @@ -208,8 +209,8 @@ fvec = fvec(:); nfev ++; else - fjac = __fdjac__ (fcn, reshape (x, xsiz), fvec); - nfev += length (x); + fjac = __fdjac__ (fcn, reshape (x, xsiz), fvec, cdif); + nfev += (1 + cdif) * length (x); endif ## For square and overdetermined systems, we update a QR