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 %%MDSMAX Multidirectional search method for direct search optimization.
18 %% [x, fmax, nf] = MDSMAX(FUN, x0, STOPIT, SAVIT) attempts to
19 %% maximize the function FUN, using the starting vector x0.
20 %% The method of multidirectional search 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 %% If a non-empty fourth parameter string SAVIT is present, then
37 %% `SAVE SAVIT x fmax nf' is executed after each inner iteration.
38 %% NB: x0 can be a matrix. In the output argument, in SAVIT saves,
39 %% and in function calls, x has the same shape as x0.
40 %% MDSMAX(fun, x0, STOPIT, SAVIT, P1, P2,...) allows additional
41 %% arguments to be passed to fun, via feval(fun,x,P1,P2,...).
43 %% This implementation uses 2n^2 elements of storage (two simplices), where x0
44 %% is an n-vector. It is based on the algorithm statement in [2, sec.3],
45 %% modified so as to halve the storage (with a slight loss in readability).
48 %% [1] V. J. Torczon, Multi-directional search: A direct search algorithm for
49 %% parallel machines, Ph.D. Thesis, Rice University, Houston, Texas, 1989.
50 % [2] V. J. Torczon, On the convergence of the multidirectional search
51 %% algorithm, SIAM J. Optimization, 1 (1991), pp. 123-145.
52 %% [3] N. J. Higham, Optimization by direct search in matrix computations,
53 %% SIAM J. Matrix Anal. Appl, 14(2): 317-333, 1993.
54 %% [4] N. J. Higham, Accuracy and Stability of Numerical Algorithms,
55 %% Second edition, Society for Industrial and Applied Mathematics,
56 %% Philadelphia, PA, 2002; sec. 20.5.
59 % Copyright (C) 2002 N.J.Higham
60 % www.maths.man.ac.uk/~higham/mctoolbox
61 % Modifications for octave by A.Adler 2003
63 function [x, fmax, nf] = mdsmax(fun, x, stopit, savit, varargin)
65 x0 = x(:); % Work with column vector internally.
68 mu = 2; % Expansion factor.
69 theta = 0.5; % Contraction factor.
71 % Set up convergence parameters etc.
74 elseif isempty(stopit)
77 tol = stopit(1); % Tolerance for cgce test based on relative size of simplex.
78 if length(stopit) == 1, stopit(2) = inf; end % Max no. of f-evaluations.
79 if length(stopit) == 2, stopit(3) = inf; end % Default target for f-values.
80 if length(stopit) == 3, stopit(4) = 0; end % Default initial simplex.
81 if length(stopit) == 4, stopit(5) = 1; end % Default: show progress.
83 if length(stopit) == 5, stopit(6) = 1; end % Default: maximize
85 if nargin < 4, savit = []; end % File name for snapshots.
87 V = [zeros(n,1) eye(n)]; T = V;
88 f = zeros(n+1,1); ft = f;
89 V(:,1) = x0; f(1) = dirn*feval(fun,x,varargin{:});
92 if trace, fprintf('f(x0) = %9.4e\n', f(1)), end
96 % Set up initial simplex.
97 scale = max(norm(x0,inf),1);
99 % Regular simplex - all edges have same length.
100 % Generated from construction given in reference [18, pp. 80-81] of [1].
101 alpha = scale / (n*sqrt(2)) * [ sqrt(n+1)-1+n sqrt(n+1)-1 ];
102 V(:,2:n+1) = (x0 + alpha(2)*ones(n,1)) * ones(1,n);
104 V(j-1,j) = x0(j-1) + alpha(1);
105 x(:) = V(:,j); f(j) = dirn*feval(fun,x,varargin{:});
108 % Right-angled simplex based on co-ordinate axes.
109 alpha = scale*ones(n+1,1);
111 V(:,j) = x0 + alpha(j)*V(:,j);
112 x(:) = V(:,j); f(j) = dirn*feval(fun,x,varargin{:});
116 size = 0; % Integer that keeps track of expansions/contractions.
117 flag_break = 0; % Flag which becomes true when ready to quit outer loop.
119 while 1 %%%%%% Outer loop.
122 % Find a new best vertex x and function value fmax = f(x).
124 V(:,[1 j]) = V(:,[j 1]); v1 = V(:,1);
125 if ~isempty(savit), x(:) = v1; eval(['save ' savit ' x fmax nf']), end
128 fprintf('Iter. %2.0f, inner = %2.0f, size = %2.0f, ', k, m, size)
129 fprintf('nf = %3.0f, f = %9.4e (%2.1f%%)\n', nf, fmax, ...
130 100*(fmax-fmax_old)/(abs(fmax_old)+eps))
134 % Stopping Test 1 - f reached target value?
136 msg = ['Exceeded target...quitting\n'];
141 while 1 %%% Inner repeat loop.
144 % Stopping Test 2 - too many f-evals?
146 msg = ['Max no. of function evaluations exceeded...quitting\n'];
147 flag_break = 1; break % Quit.
150 % Stopping Test 3 - converged? This is test (4.3) in [1].
151 size_simplex = norm(V(:,2:n+1)- v1(:,ones(1,n)),1) / max(1, norm(v1,1));
152 if size_simplex <= tol
153 msg = sprintf('Simplex size %9.4e <= %9.4e...quitting\n', ...
155 flag_break = 1; break % Quit.
158 for j=2:n+1 % ---Rotation (reflection) step.
159 T(:,j) = 2*v1 - V(:,j);
160 x(:) = T(:,j); ft(j) = dirn*feval(fun,x,varargin{:});
164 replaced = ( max(ft(2:n+1)) > fmax );
167 for j=2:n+1 % ---Expansion step.
168 V(:,j) = (1-mu)*v1 + mu*T(:,j);
169 x(:) = V(:,j); f(j) = dirn*feval(fun,x,varargin{:});
172 % Accept expansion or rotation?
173 if max(ft(2:n+1)) > max(f(2:n+1))
174 V(:,2:n+1) = T(:,2:n+1); f(2:n+1) = ft(2:n+1); % Accept rotation.
176 size = size + 1; % Accept expansion (f and V already set).
179 for j=2:n+1 % ---Contraction step.
180 V(:,j) = (1+theta)*v1 - theta*T(:,j);
181 x(:) = V(:,j); f(j) = dirn*feval(fun,x,varargin{:});
184 replaced = ( max(f(2:n+1)) > fmax );
185 % Accept contraction (f and V already set).
189 if replaced, break, end
190 if (trace && rem(m, 10) == 0)
191 fprintf(' ...inner = %2.0f...\n', m);
193 end %%% Of inner repeat loop.
195 if flag_break, break, end
196 end %%%%%% Of outer loop.
199 if trace, fprintf(msg), end