octave_packages/optim-1.2.0/mdsmax.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 %%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.
  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 %%        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,...).
  42 %%
  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).
  46 %%
  47 %% References:
  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.
  57
  58 % From Matrix Toolbox
  59 % Copyright (C) 2002 N.J.Higham
  60 % www.maths.man.ac.uk/~higham/mctoolbox
  61 % Modifications for octave by A.Adler 2003
  62
  63 function [x, fmax, nf] = mdsmax(fun, x, stopit, savit, varargin)
  64
  65 x0 = x(:);  % Work with column vector internally.
  66 n = length(x0);
  67
  68 mu = 2;      % Expansion factor.
  69 theta = 0.5; % Contraction factor.
  70
  71 % Set up convergence parameters etc.
  72 if nargin < 3
  73         stopit(1) = 1e-3;
  74 elseif isempty(stopit)
  75         stopit(1) = 1e-3;
  76 endif
  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.
  82 trace  = stopit(5);
  83 if length(stopit) == 5, stopit(6) = 1; end    % Default: maximize
  84 dirn= stopit(6);
  85 if nargin < 4, savit = []; end                   % File name for snapshots.
  86
  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{:});
  90 fmax_old = f(1);
  91
  92 if trace, fprintf('f(x0) = %9.4e\n', f(1)), end
  93
  94 k = 0; m = 0;
  95
  96 % Set up initial simplex.
  97 scale = max(norm(x0,inf),1);
  98 if stopit(4) == 0
  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);
 103    for j=2:n+1
 104        V(j-1,j) = x0(j-1) + alpha(1);
 105        x(:) = V(:,j); f(j) = dirn*feval(fun,x,varargin{:});
 106    end
 107 else
 108    % Right-angled simplex based on co-ordinate axes.
 109    alpha = scale*ones(n+1,1);
 110    for j=2:n+1
 111        V(:,j) = x0 + alpha(j)*V(:,j);
 112        x(:) = V(:,j); f(j) = dirn*feval(fun,x,varargin{:});
 113    end
 114 end
 115 nf = n+1;
 116 size = 0;         % Integer that keeps track of expansions/contractions.
 117 flag_break = 0;   % Flag which becomes true when ready to quit outer loop.
 118
 119 while 1    %%%%%% Outer loop.
 120 k = k+1;
 121
 122 % Find a new best vertex  x  and function value  fmax = f(x).
 123 [fmax,j] = max(f);
 124 V(:,[1 j]) = V(:,[j 1]); v1 = V(:,1);
 125 if ~isempty(savit), x(:) = v1; eval(['save ' savit ' x fmax nf']), end
 126 f([1 j]) = f([j 1]);
 127 if trace
 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))
 131 end
 132 fmax_old = fmax;
 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 m = 0;
 141 while 1   %%% Inner repeat loop.
 142     m = m+1;
 143
 144     % Stopping Test 2 - too many f-evals?
 145     if nf >= stopit(2)
 146        msg = ['Max no. of function evaluations exceeded...quitting\n'];
 147        flag_break = 1; break  % Quit.
 148     end
 149
 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', ...
 154                       size_simplex, tol);
 155        flag_break = 1; break  % Quit.
 156     end
 157
 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{:});
 161     end
 162     nf = nf + n;
 163
 164     replaced = ( max(ft(2:n+1)) > fmax );
 165
 166     if replaced
 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{:});
 170        end
 171        nf = nf + n;
 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.
 175        else
 176           size = size + 1;  % Accept expansion (f and V already set).
 177        end
 178     else
 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{:});
 182        end
 183        nf = nf + n;
 184        replaced = ( max(f(2:n+1)) > fmax );
 185        % Accept contraction (f and V already set).
 186        size = size - 1;
 187     end
 188
 189     if replaced, break, end
 190     if (trace && rem(m, 10) == 0)
 191       fprintf('        ...inner = %2.0f...\n', m);
 192     end
 193     end %%% Of inner repeat loop.
 194
 195 if flag_break, break, end
 196 end %%%%%% Of outer loop.
 197
 198 % Finished.
 199 if trace, fprintf(msg), end
 200 x(:) = v1;