Mercurial > octave
view scripts/optimization/fminsearch.m @ 28223:45763d59cb4f stable
use nullptr instead of NULL or 0 in a few more places
* QWinTerminalImpl.cpp, oct-procbuf.cc, audioread.cc, jit-typeinfo.cc,
lo-sysdep.cc, url-transfer.cc, shared-fcns.h: Replace NULL and 0 with
nullptr where appropriate.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Wed, 15 Apr 2020 15:55:32 -0400 |
parents | bd51beb6205e |
children | 90fea9cc9caa 0a5b15007766 |
line wrap: on
line source
######################################################################## ## ## Copyright (C) 2003-2020 The Octave Project Developers ## ## See the file COPYRIGHT.md in the top-level directory of this ## distribution or <https://octave.org/copyright/>. ## ## 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 ## <https://www.gnu.org/licenses/>. ## ######################################################################## ## -*- texinfo -*- ## @deftypefn {} {@var{x} =} fminsearch (@var{fun}, @var{x0}) ## @deftypefnx {} {@var{x} =} fminsearch (@var{fun}, @var{x0}, @var{options}) ## @deftypefnx {} {@var{x} =} fminsearch (@var{problem}) ## @deftypefnx {} {[@var{x}, @var{fval}, @var{exitflag}, @var{output}] =} fminsearch (@dots{}) ## ## Find a value of @var{x} which minimizes the multi-variable function ## @var{fun}. ## ## @var{fun} is a function handle, inline function, or string containing the ## name of the function to evaluate. ## ## The search begins at the point @var{x0} and iterates using the ## @nospell{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: ## @qcode{"Display"}, @qcode{"FunValCheck"},@qcode{"MaxFunEvals"}, ## @qcode{"MaxIter"}, @qcode{"OutputFcn"}, @qcode{"TolFun"}, @qcode{"TolX"}. ## ## @qcode{"MaxFunEvals"} proscribes the maximum number of function evaluations ## before optimization is halted. The default value is ## @code{200 * number_of_variables}, i.e., @code{200 * length (@var{x0})}. ## The value must be a positive integer. ## ## @qcode{"MaxIter"} proscribes the maximum number of algorithm iterations ## before optimization is halted. The default value is ## @code{200 * number_of_variables}, i.e., @code{200 * length (@var{x0})}. ## The value must be a positive integer. ## ## For a description of the other options, see @code{optimset}. To initialize ## an options structure with default values for @code{fminsearch} use ## @code{options = optimset ("fminsearch")}. ## ## @code{fminsearch} may also be called with a single structure argument ## with the following fields: ## ## @table @code ## @item objective ## The objective function. ## ## @item x0 ## The initial point. ## ## @item solver ## Must be set to @qcode{"fminsearch"}. ## ## @item options ## A structure returned from @code{optimset} or an empty matrix to ## indicate that defaults should be used. ## @end table ## ## @noindent ## The field @code{options} is optional. All others are required. ## ## On exit, the function returns @var{x}, the minimum point, and @var{fval}, ## the function value at the minimum. ## ## The third output @var{exitflag} reports whether the algorithm succeeded and ## may take one of the following values: ## ## @table @asis ## @item 1 ## if the algorithm converged ## (size of the simplex is smaller than @code{TolX} @strong{AND} the step in ## function value between iterations is smaller than @code{TolFun}). ## ## @item 0 ## if the maximum number of iterations or the maximum number of function ## evaluations are exceeded. ## ## @item -1 ## if the iteration is stopped by the @qcode{"OutputFcn"}. ## @end table ## ## The fourth output is a structure @var{output} containing runtime ## about the algorithm. Fields in the structure are @code{funcCount} ## containing the number of function calls to @var{fun}, @code{iterations} ## containing the number of iteration steps, @code{algorithm} with the name of ## the search algorithm (always: ## @nospell{@qcode{"Nelder-Mead simplex direct search"}}), and @code{message} ## with the exit message. ## ## Example: ## ## @example ## fminsearch (@@(x) (x(1)-5).^2+(x(2)-8).^4, [0;0]) ## @end example ## ## Note: If you need to find the minimum of a single variable function it is ## probably better to use @code{fminbnd}. ## @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 output function with "state" set to "interrupt". function [x, fval, exitflag, output] = fminsearch (varargin) if (nargin < 1) print_usage (); endif ## Get default options if requested. if (nargin == 1 && ischar (varargin{1}) && strcmp (varargin{1}, "defaults")) x = struct ("Display", "notify", "FunValCheck", "off", "MaxFunEvals", [], "MaxIter", [], "OutputFcn", [], "TolFun", 1e-4, "TolX", 1e-4); return; endif if (nargin == 1) problem = varargin{1}; varargin = {}; if (! isstruct (problem)) error ("fminsearch: PROBLEM must be a structure"); endif fun = problem.objective; x0 = problem.x0; if (! strcmp (problem.solver, "fminsearch")) error ('fminsearch: problem.solver must be set to "fminsearch"'); endif if (isfield (problem, "options")) options = problem.options; else options = []; endif elseif (nargin > 1) fun = varargin{1}; x0 = varargin{2}; if (nargin > 2) options = varargin{3}; varargin(1:3) = []; else options = []; varargin = {}; endif endif if (ischar (fun)) fun = str2func (fun); endif if (isempty (options)) options = struct (); endif [x, exitflag, output] = nmsmax (fun, x0, options, varargin{:}); 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 (i.e., 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, 2013 N.J.Higham ## www.maths.man.ac.uk/~higham/mctoolbox ## ## Modifications for Octave by A.Adler 2003 function [stopit, savit, dirn, trace, tol, maxiter, tol_f, outfcn] = ... parse_options (options, x) ## Tolerance for cgce test based on relative size of simplex. stopit(1) = tol = optimget (options, "TolX", 1e-4); ## Tolerance for cgce test based on step in function value. tol_f = optimget (options, "TolFun", 1e-4); ## Max number of function evaluations. stopit(2) = optimget (options, "MaxFunEvals", 200 * length (x)); ## Max number of iterations maxiter = optimget (options, "MaxIter", 200 * length (x)); ## Default target for function 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"); switch (display) case "iter" stopit(5) = 1; case "final" stopit(5) = 2; case "notify" stopit(5) = 3; otherwise # "none" stopit(5) = 0; endswitch trace = stopit(5); ## Use function to minimize, not maximize stopit(6) = dirn = -1; ## Filename for snapshots. savit = []; # FIXME: expose this parameter to the outside ## OutputFcn outfcn = optimget (options, "OutputFcn"); endfunction function [x, exitflag, output] = nmsmax (fun, x, options, varargin) [stopit, savit, dirn, trace, tol, maxiter, tol_f, outfcn] = ... parse_options (options, x); if (strcmpi (optimget (options, "FunValCheck", "off"), "on")) ## Replace fcn with a guarded version. fun = @(x) guarded_eval (fun, x, varargin{:}); else fun = @(x) real (fun (x, varargin{:})); 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 * fun (x); fmax_old = f(1); nf = 1; if (trace == 1) printf ("f(x0) = %9.4e\n", dirn * 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 * fun (x); 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 * fun (x); endfor endif nf += n; how = "initial "; [~, j] = sort (f); j = j(n+1:-1:1); f = f(j); V = V(:,j); exitflag = 0; if (! isempty (outfcn)) optimvalues.iteration = 0; optimvalues.funccount = nf; optimvalues.fval = dirn * f(1); optimvalues.procedure = how; state = "init"; stop = outfcn (x, optimvalues, state); if (stop) msg = "Stopped by OutputFcn\n"; exitflag = -1; endif endif alpha = 1; beta = 1/2; gamma = 2; while (exitflag != -1) # Outer (and only) loop. k += 1; if (k > maxiter) msg = "Exceeded maximum iterations\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 == 1) printf ("Iter. %2.0f,", k); printf (" how = %-11s", [how ","]); printf ("nf = %3.0f, f = %9.4e (%2.1f%%)\n", nf, dirn * 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"; ## FIXME: Add documentation when stopit(3) gets exposed to the outside exitflag = -1; break; endif ## Stopping Test 2 - too many f-evals? if (nf >= stopit(2)) msg = "Exceeded maximum number of function evaluations\n"; exitflag = 0; break; endif ## Stopping Test 3 - converged? The first part 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)); step_f = max (abs (f(1) - f(2:n+1))); if (size_simplex <= tol && step_f <= tol_f ) msg = sprintf (["Algorithm converged. Simplex size %9.4e <= %9.4e ", ... "and step in function value %9.4e <= %9.4e\n"], ... size_simplex, tol, step_f, tol_f); exitflag = 1; break; endif ## Call OutputFcn if (! isempty (outfcn)) optimvalues.funccount = nf; optimvalues.fval = dirn * f(1); optimvalues.iteration = k; optimvalues.procedure = how; state = "iter"; stop = outfcn (x, optimvalues, state); if (stop) msg = "Stopped by OutputFcn\n"; exitflag = -1; break; endif 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 * fun (x); 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 * fun (x); 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 * fun (x); 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 * fun (x); endfor nf += n-1; vk = (V(:,1) + V(:,n+1))/2; x(:) = vk; fk = dirn * fun (x); 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 == 1) || (trace == 2) || (trace == 3 && exitflag != 1) ) printf (msg); endif x(:) = V(:,1); ## FIXME: Should outputfcn be called only if exitflag != 0, ## i.e., only when we have successfully converged? if (! isempty (outfcn)) optimvalues.funccount = nf; optimvalues.fval = dirn * f(1); optimvalues.iteration = k; optimvalues.procedure = how; state = "done"; outfcn (x, optimvalues, state); endif ## output output.iterations = k; output.funcCount = nf; output.algorithm = "Nelder-Mead simplex direct search"; output.message = msg; endfunction ## A helper function that evaluates a function and checks for bad results. function y = guarded_eval (fun, x, varargin) y = fun (x, varargin{:}); if (! (isreal (y))) error ("fminsearch:notreal", "fminsearch: non-real value encountered"); elseif (any (isnan (y(:)))) error ("fminsearch:isnan", "fminsearch: NaN value encountered"); elseif (any (isinf (y(:)))) error ("fminsearch:isinf", "fminsearch: Inf value encountered"); endif endfunction %!demo %! clf; %! hold on; %! draw_fcn = @(x) (plot (x(1), x(2)) && false); %! fcn = @(x) (x(1)-5).^2 + (x(2)-8).^4; %! x0 = [0;0]; %! [xmin, fval] = fminsearch (fcn, x0, optimset ("OutputFcn", draw_fcn)) %! hold off; %!assert (fminsearch (@sin, 3, optimset ("MaxIter", 30)), 3*pi/2, 1e-4) ## FIXME: The following test is for checking that fminsearch stops earlier ## with these settings. If the optimizer algorithm is changed, it ## may fail. Just adapt the values to make it pass again. %!test %! x = fminsearch (@sin, 3, optimset ("MaxIter", 3, "Display", "none")); %! assert (x, 4.8750, 1e-4); %! x = fminsearch (@sin, 3, optimset ("MaxFunEvals", 18, "Display", "none")); %! assert (x, 4.7109, 1e-4); %!test %! problem.objective = @sin; %! problem.x0 = 3; %! problem.solver = "fminsearch"; %! problem.options = optimset ("MaxIter", 3, "Display", "none"); %! x = fminsearch (problem); %! assert (x, 4.8750, 1e-4); %! problem.options = optimset ("MaxFunEvals", 18, "Display", "none"); %! x = fminsearch (problem); %! assert (x, 4.7109, 1e-4); %!test %! c = 1.5; %! assert (fminsearch (@(x) x(1).^2 + c*x(2).^2, [1;1]), [0;0], 1e-4); ## additional input argument %!test %! x1 = fminsearch (@(x, c) x(1).^2 + c*x(2).^2, [1;1], [], 1.5); %! assert (x1, [0;0], 1e-4); %! x1 = fminsearch (@(x, c) c(1)*x(1).^2 + c(2)*x(2).^2, [1;1], ... %! optimset ("Display", "none"), [1 1.5]); %! assert (x1, [0;0], 1e-4); ## all output arguments %!test %! options = optimset ("Display", "none", "TolX", 1e-4, "TolFun", 1e-7); %! [x, fval, exitflag, output] = fminsearch (@sin, 3, options); %! assert (x, 3*pi/2, options.TolX); %! assert (fval, -1, options.TolFun); %! assert (exitflag, 1); %! assert (isstruct (output)); %! assert (isfield (output, "iterations") && isnumeric (output.iterations) %! && isscalar (output.iterations) && output.iterations > 0); %! assert (isfield (output, "funcCount") && isnumeric (output.funcCount) %! && isscalar (output.funcCount) && output.funcCount > 0); %! assert (isfield (output, "algorithm") && ischar (output.algorithm)); %! assert (isfield (output, "message") && ischar (output.message)); ## Tests for guarded_eval %!error <non-real value encountered> %! fminsearch (@(x) ([0 2i]), 0, optimset ("FunValCheck", "on")); %!error <NaN value encountered> %! fminsearch (@(x) (NaN), 0, optimset ("FunValCheck", "on")); %!error <Inf value encountered> %! fminsearch (@(x) (Inf), 0, optimset ("FunValCheck", "on")); ## Test input validation %!error fminsearch () %!error fminsearch (1)