Mercurial > octave-nkf

comparison scripts/optimization/fminsearch.m @ 14895:e0525ecf156e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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
children 6437fa7263dd
comparison
equal deleted inserted replaced
14893:55d0f8d70fe9 14895:e0525ecf156e
1 ## Copyright (C) 2003,2012 Andy Adler
2 ## Copyright (C) 2002 N.J.Higham
3 ##
4 ## This file is part of Octave.
5 ##
6 ## Octave is free software; you can redistribute it and/or modify it
7 ## under the terms of the GNU General Public License as published by
8 ## the Free Software Foundation; either version 3 of the License, or (at
9 ## your option) any later version.
10 ##
11 ## Octave is distributed in the hope that it will be useful, but
12 ## WITHOUT ANY WARRANTY; without even the implied warranty of
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 ## General Public License for more details.
15 ##
16 ## You should have received a copy of the GNU General Public License
17 ## along with Octave; see the file COPYING. If not, see
18 ## <http://www.gnu.org/licenses/>.
19
20 ## -*- texinfo -*-
21 ## @deftypefn {Function File} {@var{x} =} fminsearch (@var{fun}, @var{x0})
22 ## @deftypefnx {Function File} {@var{x} =} fminsearch (@var{fun}, @var{x0}, @var{options})
23 ## @deftypefnx {Function File} {[@var{x}, @var{fval}] =} fminsearch (@dots{})
24 ##
25 ## Find a value of @var{x} which minimizes the function @var{fun}.
26 ## The search begins at the point @var{x0} and iterates using the
27 ## Nelder & Mead Simplex algorithm (a derivative-free method). This algorithm
28 ## is better-suited to functions which have discontinuities or for which
29 ## a gradient-based search such as @code{fminunc} fails.
30 ##
31 ## Options for the search are provided in the parameter @var{options} using
32 ## the function @code{optimset}. Currently, @code{fminsearch} accepts the
33 ## options: "TolX", "MaxFunEvals", "MaxIter", "Display". For a description of
34 ## these options, see @code{optimset}.
35 ##
36 ## On exit, the function returns @var{x}, the minimum point,
37 ## and @var{fval}, the function value thereof.
38 ##
39 ## Example usages:
40 ##
41 ## @example
42 ## @group
43 ## fminsearch (@@(x) (x(1)-5).^2+(x(2)-8).^4, [0;0])
44 ##
45 ## fminsearch (inline ("(x(1)-5).^2+(x(2)-8).^4", "x"), [0;0])
46 ## @end group
47 ## @end example
48 ## @seealso{fminbnd, fminunc, optimset}
49 ## @end deftypefn
50
51 ## PKG_ADD: ## Discard result to avoid polluting workspace with ans at startup.
52 ## PKG_ADD: [~] = __all_opts__ ("fminsearch");
53
54 ## FIXME: Add support for "exitflag" output variable
55 ## FIXME: Add support for "output" output variable
56 ## FIXME: For Display option, add 'final' and 'notify' options. Not too hard.
57 ## FIXME: Add support for OutputFcn. See fminunc for a template
58 ## FIXME: Add support for exiting based on TolFun. See fminunc for an idea.
59
60 function [x, fval] = fminsearch (fun, x0, options = struct ())
61
62 ## Get default options if requested.
63 if (nargin == 1 && ischar (fun) && strcmp (fun, "defaults"))
64 x = optimset ("Display", "notify", "FunValCheck", "off",
65 "MaxFunEvals", 400, "MaxIter", 400,
66 "OutputFcn", [],
67 "TolFun", 1e-7, "TolX", 1e-4);
68 return;
69 endif
70
71 if (nargin < 2 || nargin > 3)
72 print_usage ();
73 endif
74
75 x = nmsmax (fun, x0, options);
76
77 if (isargout (2))
78 fval = feval (fun, x);
79 endif
80
81 endfunction
82
83 ##NMSMAX Nelder-Mead simplex method for direct search optimization.
84 ## [x, fmax, nf] = NMSMAX(FUN, x0, STOPIT, SAVIT) attempts to
85 ## maximize the function FUN, using the starting vector x0.
86 ## The Nelder-Mead direct search method is used.
87 ## Output arguments:
88 ## x = vector yielding largest function value found,
89 ## fmax = function value at x,
90 ## nf = number of function evaluations.
91 ## The iteration is terminated when either
92 ## - the relative size of the simplex is <= STOPIT(1)
93 ## (default 1e-3),
94 ## - STOPIT(2) function evaluations have been performed
95 ## (default inf, i.e., no limit), or
96 ## - a function value equals or exceeds STOPIT(3)
97 ## (default inf, i.e., no test on function values).
98 ## The form of the initial simplex is determined by STOPIT(4):
99 ## STOPIT(4) = 0: regular simplex (sides of equal length, the default)
100 ## STOPIT(4) = 1: right-angled simplex.
101 ## Progress of the iteration is not shown if STOPIT(5) = 0 (default 1).
102 ## STOPIT(6) indicates the direction (ie. minimization or
103 ## maximization.) Default is 1, maximization.
104 ## set STOPIT(6)=-1 for minimization
105 ## If a non-empty fourth parameter string SAVIT is present, then
106 ## `SAVE SAVIT x fmax nf' is executed after each inner iteration.
107 ## NB: x0 can be a matrix. In the output argument, in SAVIT saves,
108 ## and in function calls, x has the same shape as x0.
109 ## NMSMAX(fun, x0, STOPIT, SAVIT, P1, P2,...) allows additional
110 ## arguments to be passed to fun, via feval(fun,x,P1,P2,...).
111 ## References:
112 ## N. J. Higham, Optimization by direct search in matrix computations,
113 ## SIAM J. Matrix Anal. Appl, 14(2): 317-333, 1993.
114 ## C. T. Kelley, Iterative Methods for Optimization, Society for Industrial
115 ## and Applied Mathematics, Philadelphia, PA, 1999.
116
117 ## From Matrix Toolbox
118 ## Copyright (C) 2002 N.J.Higham
119 ## www.maths.man.ac.uk/~higham/mctoolbox
120 ##
121 ## Modifications for Octave by A.Adler 2003
122
123 function [stopit, savit, dirn, trace, tol, maxiter] = parse_options (options, x );
124
125 ## Tolerance for cgce test based on relative size of simplex.
126 stopit(1) = tol = optimget (options, "TolX", 1e-4);
127
128 ## Max no. of f-evaluations.
129 stopit(2) = optimget (options, "MaxFunEvals", length (x) * 200);
130
131 ## Max no. of iterations
132 maxiter = optimget (options, "MaxIter", length (x) * 200);
133
134 ## Default target for f-values.
135 stopit(3) = Inf; # FIXME: expose this parameter to the outside
136
137 ## Default initial simplex.
138 stopit(4) = 0; # FIXME: expose this parameter to the outside
139
140 ## Default: show progress.
141 display = optimget (options, "Display", "notify");
142 if (strcmp (display, "iter"))
143 stopit(5) = 1;
144 else
145 stopit(5) = 0;
146 endif
147 trace = stopit(5);
148
149 ## Use function to minimize, not maximize
150 stopit(6) = dirn = -1;
151
152 ## File name for snapshots.
153 savit = []; # FIXME: expose this parameter to the outside
154
155 endfunction
156
157 function [x, fmax, nf] = nmsmax (fun, x, options, savit, varargin)
158
159 [stopit, savit, dirn, trace, tol, maxiter] = parse_options (options, x);
160
161 if (strcmpi (optimget (options, "FunValCheck", "off"), "on"))
162 ## Replace fcn with a guarded version.
163 fun = @(x) guarded_eval (fun, x);
164 endif
165
166 x0 = x(:); # Work with column vector internally.
167 n = length (x0);
168
169 V = [zeros(n,1) eye(n)];
170 f = zeros (n+1,1);
171 V(:,1) = x0;
172 f(1) = dirn * feval (fun,x,varargin{:});
173 fmax_old = f(1);
174
175 if (trace)
176 fprintf ("f(x0) = %9.4e\n", f(1));
177 endif
178
179 k = 0; m = 0;
180
181 ## Set up initial simplex.
182 scale = max (norm (x0,Inf), 1);
183 if (stopit(4) == 0)
184 ## Regular simplex - all edges have same length.
185 ## Generated from construction given in reference [18, pp. 80-81] of [1].
186 alpha = scale / (n*sqrt (2)) * [sqrt(n+1)-1+n, sqrt(n+1)-1];
187 V(:,2:n+1) = (x0 + alpha(2)*ones (n,1)) * ones (1,n);
188 for j = 2:n+1
189 V(j-1,j) = x0(j-1) + alpha(1);
190 x(:) = V(:,j);
191 f(j) = dirn * feval (fun,x,varargin{:});
192 endfor
193 else
194 ## Right-angled simplex based on co-ordinate axes.
195 alpha = scale * ones(n+1,1);
196 for j=2:n+1
197 V(:,j) = x0 + alpha(j)*V(:,j);
198 x(:) = V(:,j);
199 f(j) = dirn * feval (fun,x,varargin{:});
200 endfor
201 endif
202 nf = n+1;
203 how = "initial ";
204
205 [~,j] = sort (f);
206 j = j(n+1:-1:1);
207 f = f(j);
208 V = V(:,j);
209
210 alpha = 1; beta = 1/2; gamma = 2;
211
212 while (1) # Outer (and only) loop.
213 k++;
214
215 if (k > maxiter)
216 msg = "Exceeded maximum iterations...quitting\n";
217 break;
218 endif
219
220 fmax = f(1);
221 if (fmax > fmax_old)
222 if (! isempty (savit))
223 x(:) = V(:,1);
224 eval (["save " savit " x fmax nf"]);
225 endif
226 endif
227 if (trace)
228 fprintf ("Iter. %2.0f,", k);
229 fprintf ([" how = " how " "]);
230 fprintf ("nf = %3.0f, f = %9.4e (%2.1f%%)\n", nf, fmax, ...
231 100*(fmax-fmax_old)/(abs(fmax_old)+eps));
232 endif
233 fmax_old = fmax;
234
235 ## Three stopping tests from MDSMAX.M
236
237 ## Stopping Test 1 - f reached target value?
238 if (fmax >= stopit(3))
239 msg = "Exceeded target...quitting\n";
240 break;
241 endif
242
243 ## Stopping Test 2 - too many f-evals?
244 if (nf >= stopit(2))
245 msg = "Max no. of function evaluations exceeded...quitting\n";
246 break;
247 endif
248
249 ## Stopping Test 3 - converged? This is test (4.3) in [1].
250 v1 = V(:,1);
251 size_simplex = norm (V(:,2:n+1)-v1(:,ones (1,n)),1) / max (1, norm (v1,1));
252 if (size_simplex <= tol)
253 msg = sprintf ("Simplex size %9.4e <= %9.4e...quitting\n", ...
254 size_simplex, tol);
255 break;
256 endif
257
258 ## One step of the Nelder-Mead simplex algorithm
259 ## NJH: Altered function calls and changed CNT to NF.
260 ## Changed each `fr < f(1)' type test to `>' for maximization
261 ## and re-ordered function values after sort.
262
263 vbar = (sum (V(:,1:n)')/n)'; # Mean value
264 vr = (1 + alpha)*vbar - alpha*V(:,n+1);
265 x(:) = vr;
266 fr = dirn * feval (fun,x,varargin{:});
267 nf = nf + 1;
268 vk = vr; fk = fr; how = "reflect, ";
269 if (fr > f(n))
270 if (fr > f(1))
271 ve = gamma*vr + (1-gamma)*vbar;
272 x(:) = ve;
273 fe = dirn * feval (fun,x,varargin{:});
274 nf = nf + 1;
275 if (fe > f(1))
276 vk = ve;
277 fk = fe;
278 how = "expand, ";
279 endif
280 endif
281 else
282 vt = V(:,n+1);
283 ft = f(n+1);
284 if (fr > ft)
285 vt = vr;
286 ft = fr;
287 endif
288 vc = beta*vt + (1-beta)*vbar;
289 x(:) = vc;
290 fc = dirn * feval (fun,x,varargin{:});
291 nf = nf + 1;
292 if (fc > f(n))
293 vk = vc; fk = fc;
294 how = "contract,";
295 else
296 for j = 2:n
297 V(:,j) = (V(:,1) + V(:,j))/2;
298 x(:) = V(:,j);
299 f(j) = dirn * feval (fun,x,varargin{:});
300 endfor
301 nf = nf + n-1;
302 vk = (V(:,1) + V(:,n+1))/2;
303 x(:) = vk;
304 fk = dirn * feval (fun,x,varargin{:});
305 nf = nf + 1;
306 how = "shrink, ";
307 endif
308 endif
309 V(:,n+1) = vk;
310 f(n+1) = fk;
311 [~,j] = sort(f);
312 j = j(n+1:-1:1);
313 f = f(j);
314 V = V(:,j);
315
316 endwhile # End of outer (and only) loop.
317
318 ## Finished.
319 if (trace)
320 fprintf (msg);
321 endif
322 x(:) = V(:,1);
323
324 endfunction
325
326 ## A helper function that evaluates a function and checks for bad results.
327 function y = guarded_eval (fun, x)
328
329 y = fun (x);
330
331 if (! (isreal (f)))
332 error ("fminsearch:notreal", "fminsearch: non-real value encountered");
333 elseif (any (isnan (f(:))))
334 error ("fminsearch:isnan", "fminsearch: NaN value encountered");
335 elseif (any (isinf (f(:))))
336 error ("fminsearch:isinf", "fminsearch: Inf value encountered");
337 endif
338
339 endfunction
340
341
342 %!demo
343 %! fcn = @(x) (x(1)-5).^2 + (x(2)-8).^4
344 %! x0 = [0;0];
345 %! [xmin, fval] = fminsearch (fcn, x0)
346
347 %!assert (fminsearch (@sin, 3, optimset ("MaxIter", 3)), 4.8750, 1e-4)
348 %!assert (fminsearch (@sin, 3, optimset ("MaxIter", 30)), 4.7124, 1e-4)
349 %!shared c
350 %! c = 1.5;
351 %!assert (fminsearch (@(x) x(1).^2+c*x(2).^2,[1;1]), [0;0], 1e-4)
352