Mercurial > octave-nkf
changeset 14895:e0525ecf156e
Add new function fminsearch.m
* fminsearch.m: new function.
* optimization/module.mk: Add fminsearch to build system.
* NEWS: Add fminsearch to list of new functions in 3.8.0.
* nonlin.txi, fminbnd.m, fminunc.m: Add fminsearch to documentation.
Update other optimization functions to reference fminsearch.
| author | Andy Adler <andy@analyti.ca> |
|---|---|
| date | Fri, 20 Jul 2012 09:25:37 -0700 |
| parents | 55d0f8d70fe9 |
| children | 0ba7be7fed1c |
| files | NEWS doc/interpreter/nonlin.txi scripts/optimization/fminbnd.m scripts/optimization/fminsearch.m scripts/optimization/fminunc.m scripts/optimization/module.mk |
| diffstat | 6 files changed, 385 insertions(+), 17 deletions(-) [+] |
line wrap: on
line diff
--- a/NEWS Fri Jul 20 11:40:50 2012 -0400 +++ b/NEWS Fri Jul 20 09:25:37 2012 -0700 @@ -74,10 +74,10 @@ ** Other new functions added in 3.8.0: - betaincinv lines tetramesh - colorcube rgbplot - erfcinv shrinkfaces - findfigs splinefit + betaincinv fminsearch splinefit + colorcube lines tetramesh + erfcinv rgbplot + findfigs shrinkfaces ** The default name of the Octave crash dump file is now called octave-workspace instead of octave-core.
--- a/doc/interpreter/nonlin.txi Fri Jul 20 11:40:50 2012 -0400 +++ b/doc/interpreter/nonlin.txi Fri Jul 20 09:25:37 2012 -0700 @@ -172,7 +172,11 @@ the zeroes where it crosses the x-axis. @code{fminbnd} is designed for the simpler, but very common, case of a univariate function where the interval to search is bounded. For unbounded minimization of a function with -potentially many variables use @code{fminunc}. @xref{Optimization}, for +potentially many variables use @code{fminunc} or @code{fminsearch}. The two +functions use different internal algorithms and some knowledge of the objective +function is required. For functions which can be differentiated, @code{fminunc} +is appropriate. For functions with discontinuities, or for which a gradient +search would fail, use @code{fminsearch}. @xref{Optimization}, for minimization with the presence of constraint functions. Note that searches can be made for maxima by simply inverting the objective function @tex @@ -186,3 +190,5 @@ @DOCSTRING(fminunc) +@DOCSTRING(fminsearch) +
--- a/scripts/optimization/fminbnd.m Fri Jul 20 11:40:50 2012 -0400 +++ b/scripts/optimization/fminbnd.m Fri Jul 20 09:25:37 2012 -0700 @@ -20,14 +20,13 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {[@var{x}, @var{fval}, @var{info}, @var{output}] =} fminbnd (@var{fun}, @var{a}, @var{b}, @var{options}) -## Find a minimum point of a univariate function. @var{fun} should be a -## function -## handle or name. @var{a}, @var{b} specify a starting interval. @var{options} -## is a -## structure specifying additional options. Currently, @code{fminbnd} -## recognizes these options: @code{"FunValCheck"}, @code{"OutputFcn"}, -## @code{"TolX"}, @code{"MaxIter"}, @code{"MaxFunEvals"}. -## For description of these options, see @ref{doc-optimset,,optimset}. +## Find a minimum point of a univariate function. +## +## @var{fun} should be a function handle or name. @var{a}, @var{b} specify a +## starting interval. @var{options} is a structure specifying additional +## options. Currently, @code{fminbnd} recognizes these options: +## "FunValCheck", "OutputFcn", "TolX", "MaxIter", "MaxFunEvals". For a +## description of these options, see @ref{doc-optimset,,optimset}. ## ## On exit, the function returns @var{x}, the approximate minimum point ## and @var{fval}, the function value thereof. @@ -43,7 +42,13 @@ ## @item -1 ## The algorithm has been terminated from user output function. ## @end itemize -## @seealso{optimset, fzero, fminunc} +## +## Notes: The search for a minimum is restricted to be in the interval +## bound by @var{a} and @var{b}. If you only have an initial point +## to begin searching from you will need to use an unconstrained +## minimization algorithm such as @code{fminunc} or @code{fminsearch}. +## @code{fminbnd} internally uses a Golden Section search strategy. +## @seealso{fzero, fminunc, fminsearch, optimset} ## @end deftypefn ## This is patterned after opt/fmin.f from Netlib, which in turn is taken from
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/optimization/fminsearch.m Fri Jul 20 09:25:37 2012 -0700 @@ -0,0 +1,352 @@ +## Copyright (C) 2003,2012 Andy Adler +## Copyright (C) 2002 N.J.Higham +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{x} =} fminsearch (@var{fun}, @var{x0}) +## @deftypefnx {Function File} {@var{x} =} fminsearch (@var{fun}, @var{x0}, @var{options}) +## @deftypefnx {Function File} {[@var{x}, @var{fval}] =} fminsearch (@dots{}) +## +## Find a value of @var{x} which minimizes the function @var{fun}. +## The search begins at the point @var{x0} and iterates using the +## Nelder & Mead Simplex algorithm (a derivative-free method). This algorithm +## is better-suited to functions which have discontinuities or for which +## a gradient-based search such as @code{fminunc} fails. +## +## Options for the search are provided in the parameter @var{options} using +## the function @code{optimset}. Currently, @code{fminsearch} accepts the +## options: "TolX", "MaxFunEvals", "MaxIter", "Display". For a description of +## these options, see @code{optimset}. +## +## On exit, the function returns @var{x}, the minimum point, +## and @var{fval}, the function value thereof. +## +## Example usages: +## +## @example +## @group +## fminsearch (@@(x) (x(1)-5).^2+(x(2)-8).^4, [0;0]) +## +## fminsearch (inline ("(x(1)-5).^2+(x(2)-8).^4", "x"), [0;0]) +## @end group +## @end example +## @seealso{fminbnd, fminunc, optimset} +## @end deftypefn + +## PKG_ADD: ## Discard result to avoid polluting workspace with ans at startup. +## PKG_ADD: [~] = __all_opts__ ("fminsearch"); + +## FIXME: Add support for "exitflag" output variable +## FIXME: Add support for "output" output variable +## FIXME: For Display option, add 'final' and 'notify' options. Not too hard. +## FIXME: Add support for OutputFcn. See fminunc for a template +## FIXME: Add support for exiting based on TolFun. See fminunc for an idea. + +function [x, fval] = fminsearch (fun, x0, options = struct ()) + + ## Get default options if requested. + if (nargin == 1 && ischar (fun) && strcmp (fun, "defaults")) + x = optimset ("Display", "notify", "FunValCheck", "off", + "MaxFunEvals", 400, "MaxIter", 400, + "OutputFcn", [], + "TolFun", 1e-7, "TolX", 1e-4); + return; + endif + + if (nargin < 2 || nargin > 3) + print_usage (); + endif + + x = nmsmax (fun, x0, options); + + if (isargout (2)) + fval = feval (fun, x); + endif + +endfunction + +##NMSMAX Nelder-Mead simplex method for direct search optimization. +## [x, fmax, nf] = NMSMAX(FUN, x0, STOPIT, SAVIT) attempts to +## maximize the function FUN, using the starting vector x0. +## The Nelder-Mead direct search method is used. +## Output arguments: +## x = vector yielding largest function value found, +## fmax = function value at x, +## nf = number of function evaluations. +## The iteration is terminated when either +## - the relative size of the simplex is <= STOPIT(1) +## (default 1e-3), +## - STOPIT(2) function evaluations have been performed +## (default inf, i.e., no limit), or +## - a function value equals or exceeds STOPIT(3) +## (default inf, i.e., no test on function values). +## The form of the initial simplex is determined by STOPIT(4): +## STOPIT(4) = 0: regular simplex (sides of equal length, the default) +## STOPIT(4) = 1: right-angled simplex. +## Progress of the iteration is not shown if STOPIT(5) = 0 (default 1). +## STOPIT(6) indicates the direction (ie. minimization or +## maximization.) Default is 1, maximization. +## set STOPIT(6)=-1 for minimization +## If a non-empty fourth parameter string SAVIT is present, then +## `SAVE SAVIT x fmax nf' is executed after each inner iteration. +## NB: x0 can be a matrix. In the output argument, in SAVIT saves, +## and in function calls, x has the same shape as x0. +## NMSMAX(fun, x0, STOPIT, SAVIT, P1, P2,...) allows additional +## arguments to be passed to fun, via feval(fun,x,P1,P2,...). +## References: +## N. J. Higham, Optimization by direct search in matrix computations, +## SIAM J. Matrix Anal. Appl, 14(2): 317-333, 1993. +## C. T. Kelley, Iterative Methods for Optimization, Society for Industrial +## and Applied Mathematics, Philadelphia, PA, 1999. + +## From Matrix Toolbox +## Copyright (C) 2002 N.J.Higham +## www.maths.man.ac.uk/~higham/mctoolbox +## +## Modifications for Octave by A.Adler 2003 + +function [stopit, savit, dirn, trace, tol, maxiter] = parse_options (options, x ); + + ## Tolerance for cgce test based on relative size of simplex. + stopit(1) = tol = optimget (options, "TolX", 1e-4); + + ## Max no. of f-evaluations. + stopit(2) = optimget (options, "MaxFunEvals", length (x) * 200); + + ## Max no. of iterations + maxiter = optimget (options, "MaxIter", length (x) * 200); + + ## Default target for f-values. + stopit(3) = Inf; # FIXME: expose this parameter to the outside + + ## Default initial simplex. + stopit(4) = 0; # FIXME: expose this parameter to the outside + + ## Default: show progress. + display = optimget (options, "Display", "notify"); + if (strcmp (display, "iter")) + stopit(5) = 1; + else + stopit(5) = 0; + endif + trace = stopit(5); + + ## Use function to minimize, not maximize + stopit(6) = dirn = -1; + + ## File name for snapshots. + savit = []; # FIXME: expose this parameter to the outside + +endfunction + +function [x, fmax, nf] = nmsmax (fun, x, options, savit, varargin) + + [stopit, savit, dirn, trace, tol, maxiter] = parse_options (options, x); + + if (strcmpi (optimget (options, "FunValCheck", "off"), "on")) + ## Replace fcn with a guarded version. + fun = @(x) guarded_eval (fun, x); + endif + + x0 = x(:); # Work with column vector internally. + n = length (x0); + + V = [zeros(n,1) eye(n)]; + f = zeros (n+1,1); + V(:,1) = x0; + f(1) = dirn * feval (fun,x,varargin{:}); + fmax_old = f(1); + + if (trace) + fprintf ("f(x0) = %9.4e\n", f(1)); + endif + + k = 0; m = 0; + + ## Set up initial simplex. + scale = max (norm (x0,Inf), 1); + if (stopit(4) == 0) + ## Regular simplex - all edges have same length. + ## Generated from construction given in reference [18, pp. 80-81] of [1]. + alpha = scale / (n*sqrt (2)) * [sqrt(n+1)-1+n, sqrt(n+1)-1]; + V(:,2:n+1) = (x0 + alpha(2)*ones (n,1)) * ones (1,n); + for j = 2:n+1 + V(j-1,j) = x0(j-1) + alpha(1); + x(:) = V(:,j); + f(j) = dirn * feval (fun,x,varargin{:}); + endfor + else + ## Right-angled simplex based on co-ordinate axes. + alpha = scale * ones(n+1,1); + for j=2:n+1 + V(:,j) = x0 + alpha(j)*V(:,j); + x(:) = V(:,j); + f(j) = dirn * feval (fun,x,varargin{:}); + endfor + endif + nf = n+1; + how = "initial "; + + [~,j] = sort (f); + j = j(n+1:-1:1); + f = f(j); + V = V(:,j); + + alpha = 1; beta = 1/2; gamma = 2; + + while (1) # Outer (and only) loop. + k++; + + if (k > maxiter) + msg = "Exceeded maximum iterations...quitting\n"; + break; + endif + + fmax = f(1); + if (fmax > fmax_old) + if (! isempty (savit)) + x(:) = V(:,1); + eval (["save " savit " x fmax nf"]); + endif + endif + if (trace) + fprintf ("Iter. %2.0f,", k); + fprintf ([" how = " how " "]); + fprintf ("nf = %3.0f, f = %9.4e (%2.1f%%)\n", nf, fmax, ... + 100*(fmax-fmax_old)/(abs(fmax_old)+eps)); + endif + fmax_old = fmax; + + ## Three stopping tests from MDSMAX.M + + ## Stopping Test 1 - f reached target value? + if (fmax >= stopit(3)) + msg = "Exceeded target...quitting\n"; + break; + endif + + ## Stopping Test 2 - too many f-evals? + if (nf >= stopit(2)) + msg = "Max no. of function evaluations exceeded...quitting\n"; + break; + endif + + ## Stopping Test 3 - converged? This is test (4.3) in [1]. + v1 = V(:,1); + size_simplex = norm (V(:,2:n+1)-v1(:,ones (1,n)),1) / max (1, norm (v1,1)); + if (size_simplex <= tol) + msg = sprintf ("Simplex size %9.4e <= %9.4e...quitting\n", ... + size_simplex, tol); + break; + endif + + ## One step of the Nelder-Mead simplex algorithm + ## NJH: Altered function calls and changed CNT to NF. + ## Changed each `fr < f(1)' type test to `>' for maximization + ## and re-ordered function values after sort. + + vbar = (sum (V(:,1:n)')/n)'; # Mean value + vr = (1 + alpha)*vbar - alpha*V(:,n+1); + x(:) = vr; + fr = dirn * feval (fun,x,varargin{:}); + nf = nf + 1; + vk = vr; fk = fr; how = "reflect, "; + if (fr > f(n)) + if (fr > f(1)) + ve = gamma*vr + (1-gamma)*vbar; + x(:) = ve; + fe = dirn * feval (fun,x,varargin{:}); + nf = nf + 1; + if (fe > f(1)) + vk = ve; + fk = fe; + how = "expand, "; + endif + endif + else + vt = V(:,n+1); + ft = f(n+1); + if (fr > ft) + vt = vr; + ft = fr; + endif + vc = beta*vt + (1-beta)*vbar; + x(:) = vc; + fc = dirn * feval (fun,x,varargin{:}); + nf = nf + 1; + if (fc > f(n)) + vk = vc; fk = fc; + how = "contract,"; + else + for j = 2:n + V(:,j) = (V(:,1) + V(:,j))/2; + x(:) = V(:,j); + f(j) = dirn * feval (fun,x,varargin{:}); + endfor + nf = nf + n-1; + vk = (V(:,1) + V(:,n+1))/2; + x(:) = vk; + fk = dirn * feval (fun,x,varargin{:}); + nf = nf + 1; + how = "shrink, "; + endif + endif + V(:,n+1) = vk; + f(n+1) = fk; + [~,j] = sort(f); + j = j(n+1:-1:1); + f = f(j); + V = V(:,j); + + endwhile # End of outer (and only) loop. + + ## Finished. + if (trace) + fprintf (msg); + endif + x(:) = V(:,1); + +endfunction + +## A helper function that evaluates a function and checks for bad results. +function y = guarded_eval (fun, x) + + y = fun (x); + + if (! (isreal (f))) + error ("fminsearch:notreal", "fminsearch: non-real value encountered"); + elseif (any (isnan (f(:)))) + error ("fminsearch:isnan", "fminsearch: NaN value encountered"); + elseif (any (isinf (f(:)))) + error ("fminsearch:isinf", "fminsearch: Inf value encountered"); + endif + +endfunction + + +%!demo +%! fcn = @(x) (x(1)-5).^2 + (x(2)-8).^4 +%! x0 = [0;0]; +%! [xmin, fval] = fminsearch (fcn, x0) + +%!assert (fminsearch (@sin, 3, optimset ("MaxIter", 3)), 4.8750, 1e-4) +%!assert (fminsearch (@sin, 3, optimset ("MaxIter", 30)), 4.7124, 1e-4) +%!shared c +%! c = 1.5; +%!assert (fminsearch (@(x) x(1).^2+c*x(2).^2,[1;1]), [0;0], 1e-4) +
--- a/scripts/optimization/fminunc.m Fri Jul 20 11:40:50 2012 -0400 +++ b/scripts/optimization/fminunc.m Fri Jul 20 09:25:37 2012 -0700 @@ -24,6 +24,7 @@ ## @deftypefnx {Function File} {[@var{x}, @var{fvec}, @var{info}, @var{output}, @var{grad}, @var{hess}] =} fminunc (@var{fcn}, @dots{}) ## Solve an unconstrained optimization problem defined by the function ## @var{fcn}. +## ## @var{fcn} should accepts a vector (array) defining the unknown variables, ## and return the objective function value, optionally with gradient. ## In other words, this function attempts to determine a vector @var{x} such @@ -72,9 +73,12 @@ ## (@var{output}), the output gradient (@var{grad}) and approximate Hessian ## (@var{hess}). ## -## Note: If you only have a single nonlinear equation of one variable, using -## @code{fminbnd} is usually a much better idea. -## @seealso{fminbnd, optimset} +## Notes: If you only have a single nonlinear equation of one variable then +## using @code{fminbnd} is usually a much better idea. The algorithm used is a +## gradient search which depends on the objective function being differentiable. +## If the function has discontinuities it may be better to use a derivative-free +## algorithm such as @code{fminsearch}. +## @seealso{fminbnd, fminsearch, optimset} ## @end deftypefn ## PKG_ADD: ## Discard result to avoid polluting workspace with ans at startup.
--- a/scripts/optimization/module.mk Fri Jul 20 11:40:50 2012 -0400 +++ b/scripts/optimization/module.mk Fri Jul 20 09:25:37 2012 -0700 @@ -6,6 +6,7 @@ optimization_FCN_FILES = \ optimization/__all_opts__.m \ optimization/fminbnd.m \ + optimization/fminsearch.m \ optimization/fminunc.m \ optimization/fsolve.m \ optimization/fzero.m \