1######################################################################## 2## 3## Copyright (C) 2003-2021 The Octave Project Developers 4## 5## See the file COPYRIGHT.md in the top-level directory of this 6## distribution or <https://octave.org/copyright/>. 7## 8## This file is part of Octave. 9## 10## Octave is free software: you can redistribute it and/or modify it 11## under the terms of the GNU General Public License as published by 12## the Free Software Foundation, either version 3 of the License, or 13## (at your option) any later version. 14## 15## Octave is distributed in the hope that it will be useful, but 16## WITHOUT ANY WARRANTY; without even the implied warranty of 17## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18## GNU General Public License for more details. 19## 20## You should have received a copy of the GNU General Public License 21## along with Octave; see the file COPYING. If not, see 22## <https://www.gnu.org/licenses/>. 23## 24######################################################################## 25 26## -*- texinfo -*- 27## @deftypefn {} {@var{x} =} fminsearch (@var{fun}, @var{x0}) 28## @deftypefnx {} {@var{x} =} fminsearch (@var{fun}, @var{x0}, @var{options}) 29## @deftypefnx {} {@var{x} =} fminsearch (@var{problem}) 30## @deftypefnx {} {[@var{x}, @var{fval}, @var{exitflag}, @var{output}] =} fminsearch (@dots{}) 31## 32## Find a value of @var{x} which minimizes the multi-variable function 33## @var{fun}. 34## 35## @var{fun} is a function handle, inline function, or string containing the 36## name of the function to evaluate. 37## 38## The search begins at the point @var{x0} and iterates using the 39## @nospell{Nelder & Mead} Simplex algorithm (a derivative-free method). This 40## algorithm is better-suited to functions which have discontinuities or for 41## which a gradient-based search such as @code{fminunc} fails. 42## 43## Options for the search are provided in the parameter @var{options} using the 44## function @code{optimset}. Currently, @code{fminsearch} accepts the options: 45## @qcode{"Display"}, @qcode{"FunValCheck"},@qcode{"MaxFunEvals"}, 46## @qcode{"MaxIter"}, @qcode{"OutputFcn"}, @qcode{"TolFun"}, @qcode{"TolX"}. 47## 48## @qcode{"MaxFunEvals"} proscribes the maximum number of function evaluations 49## before optimization is halted. The default value is 50## @code{200 * number_of_variables}, i.e., @code{200 * length (@var{x0})}. 51## The value must be a positive integer. 52## 53## @qcode{"MaxIter"} proscribes the maximum number of algorithm iterations 54## before optimization is halted. The default value is 55## @code{200 * number_of_variables}, i.e., @code{200 * length (@var{x0})}. 56## The value must be a positive integer. 57## 58## For a description of the other options, see @code{optimset}. To initialize 59## an options structure with default values for @code{fminsearch} use 60## @code{options = optimset ("fminsearch")}. 61## 62## @code{fminsearch} may also be called with a single structure argument 63## with the following fields: 64## 65## @table @code 66## @item objective 67## The objective function. 68## 69## @item x0 70## The initial point. 71## 72## @item solver 73## Must be set to @qcode{"fminsearch"}. 74## 75## @item options 76## A structure returned from @code{optimset} or an empty matrix to 77## indicate that defaults should be used. 78## @end table 79## 80## @noindent 81## The field @code{options} is optional. All others are required. 82## 83## On exit, the function returns @var{x}, the minimum point, and @var{fval}, 84## the function value at the minimum. 85## 86## The third output @var{exitflag} reports whether the algorithm succeeded and 87## may take one of the following values: 88## 89## @table @asis 90## @item 1 91## if the algorithm converged 92## (size of the simplex is smaller than @code{TolX} @strong{AND} the step in 93## function value between iterations is smaller than @code{TolFun}). 94## 95## @item 0 96## if the maximum number of iterations or the maximum number of function 97## evaluations are exceeded. 98## 99## @item -1 100## if the iteration is stopped by the @qcode{"OutputFcn"}. 101## @end table 102## 103## The fourth output is a structure @var{output} containing runtime 104## about the algorithm. Fields in the structure are @code{funcCount} 105## containing the number of function calls to @var{fun}, @code{iterations} 106## containing the number of iteration steps, @code{algorithm} with the name of 107## the search algorithm (always: 108## @nospell{@qcode{"Nelder-Mead simplex direct search"}}), and @code{message} 109## with the exit message. 110## 111## Example: 112## 113## @example 114## fminsearch (@@(x) (x(1)-5).^2+(x(2)-8).^4, [0;0]) 115## @end example 116## 117## Note: If you need to find the minimum of a single variable function it is 118## probably better to use @code{fminbnd}. 119## @seealso{fminbnd, fminunc, optimset} 120## @end deftypefn 121 122## PKG_ADD: ## Discard result to avoid polluting workspace with ans at startup. 123## PKG_ADD: [~] = __all_opts__ ("fminsearch"); 124 125## FIXME: Add support for output function with "state" set to "interrupt". 126 127function [x, fval, exitflag, output] = fminsearch (varargin) 128 129 if (nargin < 1) 130 print_usage (); 131 endif 132 133 ## Get default options if requested. 134 if (nargin == 1 && ischar (varargin{1}) && strcmp (varargin{1}, "defaults")) 135 x = struct ("Display", "notify", "FunValCheck", "off", 136 "MaxFunEvals", [], "MaxIter", [], 137 "OutputFcn", [], 138 "TolFun", 1e-4, "TolX", 1e-4); 139 return; 140 endif 141 142 if (nargin == 1) 143 problem = varargin{1}; 144 varargin = {}; 145 if (! isstruct (problem)) 146 error ("fminsearch: PROBLEM must be a structure"); 147 endif 148 fun = problem.objective; 149 x0 = problem.x0; 150 if (! strcmp (problem.solver, "fminsearch")) 151 error ('fminsearch: problem.solver must be set to "fminsearch"'); 152 endif 153 if (isfield (problem, "options")) 154 options = problem.options; 155 else 156 options = []; 157 endif 158 elseif (nargin > 1) 159 fun = varargin{1}; 160 x0 = varargin{2}; 161 if (nargin > 2) 162 options = varargin{3}; 163 varargin(1:3) = []; 164 else 165 options = []; 166 varargin = {}; 167 endif 168 endif 169 170 if (ischar (fun)) 171 fun = str2func (fun); 172 endif 173 174 if (isempty (options)) 175 options = struct (); 176 endif 177 178 [x, exitflag, output] = nmsmax (fun, x0, options, varargin{:}); 179 180 if (isargout (2)) 181 fval = feval (fun, x); 182 endif 183 184endfunction 185 186## NMSMAX Nelder-Mead simplex method for direct search optimization. 187## [x, fmax, nf] = NMSMAX(FUN, x0, STOPIT, SAVIT) attempts to 188## maximize the function FUN, using the starting vector x0. 189## The Nelder-Mead direct search method is used. 190## Output arguments: 191## x = vector yielding largest function value found, 192## fmax = function value at x, 193## nf = number of function evaluations. 194## The iteration is terminated when either 195## - the relative size of the simplex is <= STOPIT(1) 196## (default 1e-3), 197## - STOPIT(2) function evaluations have been performed 198## (default inf, i.e., no limit), or 199## - a function value equals or exceeds STOPIT(3) 200## (default inf, i.e., no test on function values). 201## The form of the initial simplex is determined by STOPIT(4): 202## STOPIT(4) = 0: regular simplex (sides of equal length, the default) 203## STOPIT(4) = 1: right-angled simplex. 204## Progress of the iteration is not shown if STOPIT(5) = 0 (default 1). 205## STOPIT(6) indicates the direction (i.e., minimization or 206## maximization.) Default is 1, maximization. 207## set STOPIT(6)=-1 for minimization 208## If a non-empty fourth parameter string SAVIT is present, then 209## 'SAVE SAVIT x fmax nf' is executed after each inner iteration. 210## NB: x0 can be a matrix. In the output argument, in SAVIT saves, 211## and in function calls, x has the same shape as x0. 212## NMSMAX(fun, x0, STOPIT, SAVIT, P1, P2,...) allows additional 213## arguments to be passed to fun, via feval(fun,x,P1,P2,...). 214## References: 215## N. J. Higham, Optimization by direct search in matrix computations, 216## SIAM J. Matrix Anal. Appl, 14(2): 317-333, 1993. 217## C. T. Kelley, Iterative Methods for Optimization, Society for Industrial 218## and Applied Mathematics, Philadelphia, PA, 1999. 219 220## From Matrix Toolbox 221## Copyright (C) 2002, 2013 N.J.Higham 222## www.maths.man.ac.uk/~higham/mctoolbox 223## 224## Modifications for Octave by A.Adler 2003 225 226function [stopit, savit, dirn, trace, tol, maxiter, tol_f, outfcn] = ... 227 parse_options (options, x) 228 229 ## Tolerance for cgce test based on relative size of simplex. 230 stopit(1) = tol = optimget (options, "TolX", 1e-4); 231 232 ## Tolerance for cgce test based on step in function value. 233 tol_f = optimget (options, "TolFun", 1e-4); 234 235 ## Max number of function evaluations. 236 stopit(2) = optimget (options, "MaxFunEvals", 200 * length (x)); 237 238 ## Max number of iterations 239 maxiter = optimget (options, "MaxIter", 200 * length (x)); 240 241 ## Default target for function values. 242 stopit(3) = Inf; # FIXME: expose this parameter to the outside 243 244 ## Default initial simplex. 245 stopit(4) = 0; # FIXME: expose this parameter to the outside 246 247 ## Default: show progress. 248 display = optimget (options, "Display", "notify"); 249 switch (display) 250 case "iter" 251 stopit(5) = 1; 252 case "final" 253 stopit(5) = 2; 254 case "notify" 255 stopit(5) = 3; 256 otherwise # "none" 257 stopit(5) = 0; 258 endswitch 259 trace = stopit(5); 260 261 ## Use function to minimize, not maximize 262 stopit(6) = dirn = -1; 263 264 ## Filename for snapshots. 265 savit = []; # FIXME: expose this parameter to the outside 266 267 ## OutputFcn 268 outfcn = optimget (options, "OutputFcn"); 269 270endfunction 271 272function [x, exitflag, output] = nmsmax (fun, x, options, varargin) 273 274 [stopit, savit, dirn, trace, tol, maxiter, tol_f, outfcn] = ... 275 parse_options (options, x); 276 277 if (strcmpi (optimget (options, "FunValCheck", "off"), "on")) 278 ## Replace fcn with a guarded version. 279 fun = @(x) guarded_eval (fun, x, varargin{:}); 280 else 281 fun = @(x) real (fun (x, varargin{:})); 282 endif 283 284 x0 = x(:); # Work with column vector internally. 285 n = length (x0); 286 287 V = [zeros(n,1) eye(n)]; 288 f = zeros (n+1,1); 289 V(:,1) = x0; 290 f(1) = dirn * fun (x); 291 fmax_old = f(1); 292 nf = 1; 293 294 if (trace == 1) 295 printf ("f(x0) = %9.4e\n", dirn * f(1)); 296 endif 297 298 k = 0; m = 0; 299 300 ## Set up initial simplex. 301 scale = max (norm (x0, Inf), 1); 302 if (stopit(4) == 0) 303 ## Regular simplex - all edges have same length. 304 ## Generated from construction given in reference [18, pp. 80-81] of [1]. 305 alpha = scale / (n*sqrt (2)) * [sqrt(n+1)-1+n, sqrt(n+1)-1]; 306 V(:,2:n+1) = (x0 + alpha(2)*ones (n,1)) * ones (1,n); 307 for j = 2:n+1 308 V(j-1,j) = x0(j-1) + alpha(1); 309 x(:) = V(:,j); 310 f(j) = dirn * fun (x); 311 endfor 312 else 313 ## Right-angled simplex based on co-ordinate axes. 314 alpha = scale * ones(n+1,1); 315 for j = 2:n+1 316 V(:,j) = x0 + alpha(j)*V(:,j); 317 x(:) = V(:,j); 318 f(j) = dirn * fun (x); 319 endfor 320 endif 321 nf += n; 322 how = "initial "; 323 324 [~, j] = sort (f); 325 j = j(n+1:-1:1); 326 f = f(j); 327 V = V(:,j); 328 329 exitflag = 0; 330 331 if (! isempty (outfcn)) 332 optimvalues.iteration = 0; 333 optimvalues.funccount = nf; 334 optimvalues.fval = dirn * f(1); 335 optimvalues.procedure = how; 336 state = "init"; 337 stop = outfcn (x, optimvalues, state); 338 if (stop) 339 msg = "Stopped by OutputFcn\n"; 340 exitflag = -1; 341 endif 342 endif 343 344 alpha = 1; beta = 1/2; gamma = 2; 345 346 while (exitflag != -1) # Outer (and only) loop. 347 k += 1; 348 349 if (k > maxiter) 350 msg = "Exceeded maximum iterations\n"; 351 break; 352 endif 353 354 fmax = f(1); 355 if (fmax > fmax_old) 356 if (! isempty (savit)) 357 x(:) = V(:,1); 358 eval (["save " savit " x fmax nf"]); 359 endif 360 endif 361 if (trace == 1) 362 printf ("Iter. %2.0f,", k); 363 printf (" how = %-11s", [how ","]); 364 printf ("nf = %3.0f, f = %9.4e (%2.1f%%)\n", nf, dirn * fmax, ... 365 100*(fmax-fmax_old)/(abs(fmax_old)+eps)); 366 endif 367 fmax_old = fmax; 368 369 ## Three stopping tests from MDSMAX.M 370 371 ## Stopping Test 1 - f reached target value? 372 if (fmax >= stopit(3)) 373 msg = "Exceeded target...quitting\n"; 374 ## FIXME: Add documentation when stopit(3) gets exposed to the outside 375 exitflag = -1; 376 break; 377 endif 378 379 ## Stopping Test 2 - too many f-evals? 380 if (nf >= stopit(2)) 381 msg = "Exceeded maximum number of function evaluations\n"; 382 exitflag = 0; 383 break; 384 endif 385 386 ## Stopping Test 3 - converged? The first part is test (4.3) in [1]. 387 v1 = V(:,1); 388 size_simplex = norm (V(:,2:n+1)-v1(:,ones (1,n)),1) / max (1, norm (v1,1)); 389 step_f = max (abs (f(1) - f(2:n+1))); 390 if (size_simplex <= tol && step_f <= tol_f ) 391 msg = sprintf (["Algorithm converged. Simplex size %9.4e <= %9.4e ", ... 392 "and step in function value %9.4e <= %9.4e\n"], ... 393 size_simplex, tol, step_f, tol_f); 394 exitflag = 1; 395 break; 396 endif 397 398 ## Call OutputFcn 399 if (! isempty (outfcn)) 400 optimvalues.funccount = nf; 401 optimvalues.fval = dirn * f(1); 402 optimvalues.iteration = k; 403 optimvalues.procedure = how; 404 state = "iter"; 405 stop = outfcn (x, optimvalues, state); 406 if (stop) 407 msg = "Stopped by OutputFcn\n"; 408 exitflag = -1; 409 break; 410 endif 411 endif 412 413 ## One step of the Nelder-Mead simplex algorithm 414 ## NJH: Altered function calls and changed CNT to NF. 415 ## Changed each 'fr < f(1)' type test to '>' for maximization 416 ## and re-ordered function values after sort. 417 418 vbar = (sum (V(:,1:n)')/n)'; # Mean value 419 vr = (1 + alpha)*vbar - alpha*V(:,n+1); 420 x(:) = vr; 421 fr = dirn * fun (x); 422 nf += 1; 423 vk = vr; fk = fr; how = "reflect"; 424 if (fr > f(n)) 425 if (fr > f(1)) 426 ve = gamma*vr + (1-gamma)*vbar; 427 x(:) = ve; 428 fe = dirn * fun (x); 429 nf += 1; 430 if (fe > f(1)) 431 vk = ve; 432 fk = fe; 433 how = "expand"; 434 endif 435 endif 436 else 437 vt = V(:,n+1); 438 ft = f(n+1); 439 if (fr > ft) 440 vt = vr; 441 ft = fr; 442 endif 443 vc = beta*vt + (1-beta)*vbar; 444 x(:) = vc; 445 fc = dirn * fun (x); 446 nf += 1; 447 if (fc > f(n)) 448 vk = vc; fk = fc; 449 how = "contract"; 450 else 451 for j = 2:n 452 V(:,j) = (V(:,1) + V(:,j))/2; 453 x(:) = V(:,j); 454 f(j) = dirn * fun (x); 455 endfor 456 nf += n-1; 457 vk = (V(:,1) + V(:,n+1))/2; 458 x(:) = vk; 459 fk = dirn * fun (x); 460 nf += 1; 461 how = "shrink"; 462 endif 463 endif 464 V(:,n+1) = vk; 465 f(n+1) = fk; 466 [~,j] = sort(f); 467 j = j(n+1:-1:1); 468 f = f(j); 469 V = V(:,j); 470 471 endwhile # End of outer (and only) loop. 472 473 ## Finished. 474 if ( (trace == 1) || (trace == 2) || (trace == 3 && exitflag != 1) ) 475 printf (msg); 476 endif 477 x(:) = V(:,1); 478 479 ## FIXME: Should outputfcn be called only if exitflag != 0, 480 ## i.e., only when we have successfully converged? 481 if (! isempty (outfcn)) 482 optimvalues.funccount = nf; 483 optimvalues.fval = dirn * f(1); 484 optimvalues.iteration = k; 485 optimvalues.procedure = how; 486 state = "done"; 487 outfcn (x, optimvalues, state); 488 endif 489 490 ## output 491 output.iterations = k; 492 output.funcCount = nf; 493 output.algorithm = "Nelder-Mead simplex direct search"; 494 output.message = msg; 495 496endfunction 497 498## A helper function that evaluates a function and checks for bad results. 499function y = guarded_eval (fun, x, varargin) 500 501 y = fun (x, varargin{:}); 502 503 if (! (isreal (y))) 504 error ("fminsearch:notreal", "fminsearch: non-real value encountered"); 505 elseif (any (isnan (y(:)))) 506 error ("fminsearch:isnan", "fminsearch: NaN value encountered"); 507 elseif (any (isinf (y(:)))) 508 error ("fminsearch:isinf", "fminsearch: Inf value encountered"); 509 endif 510 511endfunction 512 513 514%!demo 515%! clf; 516%! hold on; 517%! draw_fcn = @(x) (plot (x(1), x(2)) && false); 518%! fcn = @(x) (x(1)-5).^2 + (x(2)-8).^4; 519%! x0 = [0;0]; 520%! [xmin, fval] = fminsearch (fcn, x0, optimset ("OutputFcn", draw_fcn)) 521%! hold off; 522 523%!assert (fminsearch (@sin, 3, optimset ("MaxIter", 30)), 3*pi/2, 1e-4) 524 525## FIXME: The following test is for checking that fminsearch stops earlier 526## with these settings. If the optimizer algorithm is changed, it 527## may fail. Just adapt the values to make it pass again. 528%!test 529%! x = fminsearch (@sin, 3, optimset ("MaxIter", 3, "Display", "none")); 530%! assert (x, 4.8750, 1e-4); 531%! x = fminsearch (@sin, 3, optimset ("MaxFunEvals", 18, "Display", "none")); 532%! assert (x, 4.7109, 1e-4); 533 534%!test 535%! problem.objective = @sin; 536%! problem.x0 = 3; 537%! problem.solver = "fminsearch"; 538%! problem.options = optimset ("MaxIter", 3, "Display", "none"); 539%! x = fminsearch (problem); 540%! assert (x, 4.8750, 1e-4); 541%! problem.options = optimset ("MaxFunEvals", 18, "Display", "none"); 542%! x = fminsearch (problem); 543%! assert (x, 4.7109, 1e-4); 544 545%!test 546%! c = 1.5; 547%! assert (fminsearch (@(x) x(1).^2 + c*x(2).^2, [1;1]), [0;0], 1e-4); 548 549## additional input argument 550%!test 551%! x1 = fminsearch (@(x, c) x(1).^2 + c*x(2).^2, [1;1], [], 1.5); 552%! assert (x1, [0;0], 1e-4); 553%! x1 = fminsearch (@(x, c) c(1)*x(1).^2 + c(2)*x(2).^2, [1;1], ... 554%! optimset ("Display", "none"), [1 1.5]); 555%! assert (x1, [0;0], 1e-4); 556 557## all output arguments 558%!test 559%! options = optimset ("Display", "none", "TolX", 1e-4, "TolFun", 1e-7); 560%! [x, fval, exitflag, output] = fminsearch (@sin, 3, options); 561%! assert (x, 3*pi/2, options.TolX); 562%! assert (fval, -1, options.TolFun); 563%! assert (exitflag, 1); 564%! assert (isstruct (output)); 565%! assert (isfield (output, "iterations") && isnumeric (output.iterations) 566%! && isscalar (output.iterations) && output.iterations > 0); 567%! assert (isfield (output, "funcCount") && isnumeric (output.funcCount) 568%! && isscalar (output.funcCount) && output.funcCount > 0); 569%! assert (isfield (output, "algorithm") && ischar (output.algorithm)); 570%! assert (isfield (output, "message") && ischar (output.message)); 571 572## Tests for guarded_eval 573%!error <non-real value encountered> 574%! fminsearch (@(x) ([0 2i]), 0, optimset ("FunValCheck", "on")); 575%!error <NaN value encountered> 576%! fminsearch (@(x) (NaN), 0, optimset ("FunValCheck", "on")); 577%!error <Inf value encountered> 578%! fminsearch (@(x) (Inf), 0, optimset ("FunValCheck", "on")); 579 580## Test input validation 581%!error fminsearch () 582%!error fminsearch (1) 583