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