1 %% Copyright (C) 2002 N.J.Higham
2 %% Copyright (C) 2003 Andy Adler <adler@ncf.ca>
4 %% This program is free software; you can redistribute it and/or modify it under
5 %% the terms of the GNU General Public License as published by the Free Software
6 %% Foundation; either version 3 of the License, or (at your option) any later
9 %% This program is distributed in the hope that it will be useful, but WITHOUT
10 %% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11 %% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
14 %% You should have received a copy of the GNU General Public License along with
15 %% this program; if not, see <http://www.gnu.org/licenses/>.
17 %%NMSMAX Nelder-Mead simplex method for direct search optimization.
18 %% [x, fmax, nf] = NMSMAX(FUN, x0, STOPIT, SAVIT) attempts to
19 %% maximize the function FUN, using the starting vector x0.
20 %% The Nelder-Mead direct search method is used.
22 %% x = vector yielding largest function value found,
23 %% fmax = function value at x,
24 %% nf = number of function evaluations.
25 %% The iteration is terminated when either
26 %% - the relative size of the simplex is <= STOPIT(1)
28 %% - STOPIT(2) function evaluations have been performed
29 %% (default inf, i.e., no limit), or
30 %% - a function value equals or exceeds STOPIT(3)
31 %% (default inf, i.e., no test on function values).
32 %% The form of the initial simplex is determined by STOPIT(4):
33 %% STOPIT(4) = 0: regular simplex (sides of equal length, the default)
34 %% STOPIT(4) = 1: right-angled simplex.
35 %% Progress of the iteration is not shown if STOPIT(5) = 0 (default 1).
36 %% STOPIT(6) indicates the direction (ie. minimization or
37 %% maximization.) Default is 1, maximization.
38 %% set STOPIT(6)=-1 for minimization
39 %% If a non-empty fourth parameter string SAVIT is present, then
40 %% `SAVE SAVIT x fmax nf' is executed after each inner iteration.
41 %% NB: x0 can be a matrix. In the output argument, in SAVIT saves,
42 %% and in function calls, x has the same shape as x0.
43 %% NMSMAX(fun, x0, STOPIT, SAVIT, P1, P2,...) allows additional
44 %% arguments to be passed to fun, via feval(fun,x,P1,P2,...).
46 %% N. J. Higham, Optimization by direct search in matrix computations,
47 %% SIAM J. Matrix Anal. Appl, 14(2): 317-333, 1993.
48 %% C. T. Kelley, Iterative Methods for Optimization, Society for Industrial
49 %% and Applied Mathematics, Philadelphia, PA, 1999.
52 % Copyright (C) 2002 N.J.Higham
53 % www.maths.man.ac.uk/~higham/mctoolbox
54 % Modifications for octave by A.Adler 2003
56 function [x, fmax, nf] = nmsmax(fun, x, stopit, savit, varargin)
58 x0 = x(:); % Work with column vector internally.
61 % Set up convergence parameters etc.
62 if (nargin < 3 || isempty(stopit))
65 tol = stopit(1); % Tolerance for cgce test based on relative size of simplex.
66 if length(stopit) == 1, stopit(2) = inf; end % Max no. of f-evaluations.
67 if length(stopit) == 2, stopit(3) = inf; end % Default target for f-values.
68 if length(stopit) == 3, stopit(4) = 0; end % Default initial simplex.
69 if length(stopit) == 4, stopit(5) = 1; end % Default: show progress.
71 if length(stopit) == 5, stopit(6) = 1; end % Default: maximize
73 if nargin < 4, savit = []; end % File name for snapshots.
75 V = [zeros(n,1) eye(n)];
78 f(1) = dirn*feval(fun,x,varargin{:});
81 if trace, fprintf('f(x0) = %9.4e\n', f(1)), end
85 % Set up initial simplex.
86 scale = max(norm(x0,inf),1);
88 % Regular simplex - all edges have same length.
89 % Generated from construction given in reference [18, pp. 80-81] of [1].
90 alpha = scale / (n*sqrt(2)) * [ sqrt(n+1)-1+n sqrt(n+1)-1 ];
91 V(:,2:n+1) = (x0 + alpha(2)*ones(n,1)) * ones(1,n);
93 V(j-1,j) = x0(j-1) + alpha(1);
95 f(j) = dirn*feval(fun,x,varargin{:});
98 % Right-angled simplex based on co-ordinate axes.
99 alpha = scale*ones(n+1,1);
101 V(:,j) = x0 + alpha(j)*V(:,j);
103 f(j) = dirn*feval(fun,x,varargin{:});
111 f = f(j); V = V(:,j);
113 alpha = 1; beta = 1/2; gamma = 2;
115 while 1 %%%%%% Outer (and only) loop.
121 x(:) = V(:,1); eval(['save ' savit ' x fmax nf'])
125 fprintf('Iter. %2.0f,', k)
126 fprintf([' how = ' how ' ']);
127 fprintf('nf = %3.0f, f = %9.4e (%2.1f%%)\n', nf, fmax, ...
128 100*(fmax-fmax_old)/(abs(fmax_old)+eps))
132 %%% Three stopping tests from MDSMAX.M
134 % Stopping Test 1 - f reached target value?
136 msg = ['Exceeded target...quitting\n'];
140 % Stopping Test 2 - too many f-evals?
142 msg = ['Max no. of function evaluations exceeded...quitting\n'];
146 % Stopping Test 3 - converged? This is test (4.3) in [1].
148 size_simplex = norm(V(:,2:n+1)-v1(:,ones(1,n)),1) / max(1, norm(v1,1));
149 if size_simplex <= tol
150 msg = sprintf('Simplex size %9.4e <= %9.4e...quitting\n', ...
155 % One step of the Nelder-Mead simplex algorithm
156 % NJH: Altered function calls and changed CNT to NF.
157 % Changed each `fr < f(1)' type test to `>' for maximization
158 % and re-ordered function values after sort.
160 vbar = (sum(V(:,1:n)')/n)'; % Mean value
161 vr = (1 + alpha)*vbar - alpha*V(:,n+1);
163 fr = dirn*feval(fun,x,varargin{:});
165 vk = vr; fk = fr; how = 'reflect, ';
168 ve = gamma*vr + (1-gamma)*vbar;
170 fe = dirn*feval(fun,x,varargin{:});
178 vt = V(:,n+1); ft = f(n+1);
182 vc = beta*vt + (1-beta)*vbar;
184 fc = dirn*feval(fun,x,varargin{:});
191 V(:,j) = (V(:,1) + V(:,j))/2;
193 f(j) = dirn*feval(fun,x,varargin{:});
196 vk = (V(:,1) + V(:,n+1))/2;
198 fk = dirn*feval(fun,x,varargin{:});
207 f = f(j); V = V(:,j);
209 end %%%%%% End of outer (and only) loop.
212 if trace, fprintf(msg), end