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 %%ADSMAX Alternating directions method for direct search optimization.
18 %% [x, fmax, nf] = ADSMAX(FUN, x0, STOPIT, SAVIT, P) attempts to
19 %% maximize the function FUN, using the starting vector x0.
20 %% The alternating directions 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 increase in function value between successive
27 %% iterations is <= STOPIT(1) (default 1e-3),
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 %% Progress of the iteration is not shown if STOPIT(5) = 0 (default 1).
33 %% If a non-empty fourth parameter string SAVIT is present, then
34 %% `SAVE SAVIT x fmax nf' is executed after each inner iteration.
35 %% By default, the search directions are the co-ordinate directions.
36 %% The columns of a fifth parameter matrix P specify alternative search
37 %% directions (P = EYE is the default).
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 %% ADSMAX(fun, x0, STOPIT, SAVIT, P, P1, P2,...) allows additional
41 %% arguments to be passed to fun, via feval(fun,x,P1,P2,...).
43 %% N. J. Higham, Optimization by direct search in matrix computations,
44 %% SIAM J. Matrix Anal. Appl, 14(2): 317-333, 1993.
45 %% N. J. Higham, Accuracy and Stability of Numerical Algorithms,
46 %% Second edition, Society for Industrial and Applied Mathematics,
47 %% Philadelphia, PA, 2002; sec. 20.5.
50 % Copyright (C) 2002 N.J.Higham
51 % www.maths.man.ac.uk/~higham/mctoolbox
52 % Modifications for octave by A.Adler 2003
54 function [x, fmax, nf] = adsmax(f, x, stopit, savit, P, varargin)
56 x0 = x(:); % Work with column vector internally.
59 mu = 1e-4; % Initial percentage change in components.
60 nstep = 25; % Max number of times to double or decrease h.
62 % Set up convergence parameters.
65 elseif isempty(stopit)
68 tol = stopit(1); % Required rel. increase in function value over one iteration.
69 if length(stopit) == 1, stopit(2) = inf; end % Max no. of f-evaluations.
70 if length(stopit) == 2, stopit(3) = inf; end % Default target for f-values.
71 if length(stopit) < 5, stopit(5) = 1; end % Default: show progress.
73 if length(stopit) == 5, stopit(6) = 1; end % Default: maximize
75 if nargin < 4, savit = []; end % File name for snapshots.
77 % Matrix of search directions.
83 if ~isequal(size(P),[n n]) % Check for common error.
84 error('P must be of dimension the number of elements in x0.')
88 fmax = dirn*feval(f,x,varargin{:}); nf = 1;
89 if trace, fprintf('f(x0) = %9.4e\n', fmax), end
96 if trace, fprintf('Iter %2.0f (nf = %2.0f)\n', it, nf), end
99 for i=1:n % Loop over search directions.
104 h = sign(pi'*yi)*norm(pi.*yi)*mu; % Initial step size.
105 if h == 0, h = max(norm(yi,inf),1)*mu; end
107 x(:) = y; fnew = dirn*feval(f,x,varargin{:}); nf = nf + 1;
112 fprintf('Comp. = %2.0f, steps = %2.0f, f = %9.4e*\n', i,0,fmax)
113 fprintf('Exceeded target...quitting\n')
117 h = 2*h; lim = nstep; k = 1;
119 h = -h; lim = nstep+1; k = 0;
124 x(:) = y; fnew = dirn*feval(f,x,varargin{:}); nf = nf + 1;
125 if fnew <= fmax, break, end
126 fmax = fnew; k = k + 1;
129 fprintf('Comp. = %2.0f, steps = %2.0f, f = %9.4e*\n', i,j,fmax)
130 fprintf('Exceeded target...quitting\n')
139 if k == 0, y = yi; end
142 fprintf('Comp. = %2.0f, steps = %2.0f, f = %9.4e', i, k, fmax)
143 fprintf(' (%2.1f%%)\n', 100*(fmax-flast)/(abs(flast)+eps))
149 fprintf('Max no. of function evaluations exceeded...quitting\n')
154 if (fmax > flast && ~isempty (savit))
156 eval(['save ' savit ' x fmax nf'])
159 end % Loop over search directions.
161 if isequal(steps,zeros(n,1))
162 if trace, fprintf('Stagnated...quitting\n'), end
166 if fmax-fmax_old <= tol*abs(fmax_old)
167 if trace, fprintf('Function values ''converged''...quitting\n'), end
171 end %%%%%% Of outer loop.