octave_packages/optim-1.2.0/nmsmax.m

   1 %% Copyright (C) 2002 N.J.Higham
   2 %% Copyright (C) 2003 Andy Adler <adler@ncf.ca>
   3 %%
   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
   7 %% version.
   8 %%
   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
  12 %% details.
  13 %%
  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/>.
  16
  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.
  21 %%        Output arguments:
  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)
  27 %%                 (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 %%        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,...).
  45 %% References:
  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.
  50
  51 % From Matrix Toolbox
  52 % Copyright (C) 2002 N.J.Higham
  53 % www.maths.man.ac.uk/~higham/mctoolbox
  54 % Modifications for octave by A.Adler 2003
  55
  56 function [x, fmax, nf] = nmsmax(fun, x, stopit, savit, varargin)
  57
  58 x0 = x(:);  % Work with column vector internally.
  59 n = length(x0);
  60
  61 % Set up convergence parameters etc.
  62 if (nargin < 3 || isempty(stopit))
  63   stopit(1) = 1e-3;
  64 end
  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.
  70 trace  = stopit(5);
  71 if length(stopit) == 5, stopit(6) = 1; end    % Default: maximize
  72 dirn= stopit(6);
  73 if nargin < 4, savit = []; end                   % File name for snapshots.
  74
  75 V = [zeros(n,1) eye(n)];
  76 f = zeros(n+1,1);
  77 V(:,1) = x0;
  78 f(1) = dirn*feval(fun,x,varargin{:});
  79 fmax_old = f(1);
  80
  81 if trace, fprintf('f(x0) = %9.4e\n', f(1)), end
  82
  83 k = 0; m = 0;
  84
  85 % Set up initial simplex.
  86 scale = max(norm(x0,inf),1);
  87 if stopit(4) == 0
  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);
  92    for j=2:n+1
  93        V(j-1,j) = x0(j-1) + alpha(1);
  94        x(:) = V(:,j);
  95        f(j) = dirn*feval(fun,x,varargin{:});
  96    end
  97 else
  98    % Right-angled simplex based on co-ordinate axes.
  99    alpha = scale*ones(n+1,1);
 100    for j=2:n+1
 101        V(:,j) = x0 + alpha(j)*V(:,j);
 102        x(:) = V(:,j);
 103        f(j) = dirn*feval(fun,x,varargin{:});
 104    end
 105 end
 106 nf = n+1;
 107 how = 'initial  ';
 108
 109 [temp,j] = sort(f);
 110 j = j(n+1:-1:1);
 111 f = f(j); V = V(:,j);
 112
 113 alpha = 1;  beta = 1/2;  gamma = 2;
 114
 115 while 1    %%%%%% Outer (and only) loop.
 116 k = k+1;
 117
 118     fmax = f(1);
 119     if fmax > fmax_old
 120        if ~isempty(savit)
 121           x(:) = V(:,1); eval(['save ' savit ' x fmax nf'])
 122        end
 123     end
 124     if trace
 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))
 129     end
 130     fmax_old = fmax;
 131
 132     %%% Three stopping tests from MDSMAX.M
 133
 134     % Stopping Test 1 - f reached target value?
 135     if fmax >= stopit(3)
 136        msg = ['Exceeded target...quitting\n'];
 137        break  % Quit.
 138     end
 139
 140     % Stopping Test 2 - too many f-evals?
 141     if nf >= stopit(2)
 142        msg = ['Max no. of function evaluations exceeded...quitting\n'];
 143        break  % Quit.
 144     end
 145
 146     % Stopping Test 3 - converged?   This is test (4.3) in [1].
 147     v1 = V(:,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', ...
 151                       size_simplex, tol);
 152        break  % Quit.
 153     end
 154
 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.
 159
 160     vbar = (sum(V(:,1:n)')/n)';  % Mean value
 161     vr = (1 + alpha)*vbar - alpha*V(:,n+1);
 162     x(:) = vr;
 163     fr = dirn*feval(fun,x,varargin{:});
 164     nf = nf + 1;
 165     vk = vr;  fk = fr; how = 'reflect, ';
 166     if fr > f(n)
 167         if fr > f(1)
 168            ve = gamma*vr + (1-gamma)*vbar;
 169            x(:) = ve;
 170            fe = dirn*feval(fun,x,varargin{:});
 171            nf = nf + 1;
 172            if fe > f(1)
 173               vk = ve; fk = fe;
 174               how = 'expand,  ';
 175            end
 176         end
 177     else
 178         vt = V(:,n+1); ft = f(n+1);
 179         if fr > ft
 180            vt = vr;  ft = fr;
 181         end
 182         vc = beta*vt + (1-beta)*vbar;
 183         x(:) = vc;
 184         fc = dirn*feval(fun,x,varargin{:});
 185         nf = nf + 1;
 186         if fc > f(n)
 187            vk = vc; fk = fc;
 188            how = 'contract,';
 189         else
 190            for j = 2:n
 191                V(:,j) = (V(:,1) + V(:,j))/2;
 192                x(:) = V(:,j);
 193                f(j) = dirn*feval(fun,x,varargin{:});
 194            end
 195            nf = nf + n-1;
 196            vk = (V(:,1) + V(:,n+1))/2;
 197            x(:) = vk;
 198            fk = dirn*feval(fun,x,varargin{:});
 199            nf = nf + 1;
 200            how = 'shrink,  ';
 201         end
 202     end
 203     V(:,n+1) = vk;
 204     f(n+1) = fk;
 205     [temp,j] = sort(f);
 206     j = j(n+1:-1:1);
 207     f = f(j); V = V(:,j);
 208
 209 end   %%%%%% End of outer (and only) loop.
 210
 211 % Finished.
 212 if trace, fprintf(msg), end
 213 x(:) = V(:,1);