octave_packages/optim-1.2.0/adsmax.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 %%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.
  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 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,...).
  42 %%     Reference:
  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.
  48
  49 % From Matrix Toolbox
  50 % Copyright (C) 2002 N.J.Higham
  51 % www.maths.man.ac.uk/~higham/mctoolbox
  52 % Modifications for octave by A.Adler 2003
  53
  54 function [x, fmax, nf] = adsmax(f, x, stopit, savit, P, varargin)
  55
  56 x0 = x(:);  % Work with column vector internally.
  57 n = length(x0);
  58
  59 mu = 1e-4;  % Initial percentage change in components.
  60 nstep = 25; % Max number of times to double or decrease h.
  61
  62 % Set up convergence parameters.
  63 if nargin < 3
  64         stopit(1) = 1e-3;
  65 elseif isempty(stopit)
  66         stopit(1) = 1e-3;
  67 endif
  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.
  72 trace  = stopit(5);
  73 if length(stopit) == 5, stopit(6) = 1; end    % Default: maximize
  74 dirn= stopit(6);
  75 if nargin < 4, savit = []; end                   % File name for snapshots.
  76
  77 % Matrix of search directions.
  78 if nargin < 5
  79    P = eye(n);
  80 elseif isempty(P)
  81    P = eye(n);
  82 else
  83    if ~isequal(size(P),[n n])  % Check for common error.
  84       error('P must be of dimension the number of elements in x0.')
  85    end
  86 end
  87
  88 fmax = dirn*feval(f,x,varargin{:}); nf = 1;
  89 if trace, fprintf('f(x0) = %9.4e\n', fmax), end
  90
  91 steps = zeros(n,1);
  92 it = 0; y = x0;
  93
  94 while 1    % Outer loop.
  95 it = it+1;
  96 if trace, fprintf('Iter %2.0f  (nf = %2.0f)\n', it, nf), end
  97 fmax_old = fmax;
  98
  99 for i=1:n  % Loop over search directions.
 100
 101     pi = P(:,i);
 102     flast = fmax;
 103     yi = y;
 104     h = sign(pi'*yi)*norm(pi.*yi)*mu;   % Initial step size.
 105     if h == 0, h = max(norm(yi,inf),1)*mu; end
 106     y = yi + h*pi;
 107     x(:) = y; fnew = dirn*feval(f,x,varargin{:}); nf = nf + 1;
 108     if fnew > fmax
 109        fmax = fnew;
 110        if fmax >= stopit(3)
 111            if trace
 112               fprintf('Comp. = %2.0f,  steps = %2.0f,  f = %9.4e*\n', i,0,fmax)
 113               fprintf('Exceeded target...quitting\n')
 114            end
 115            x(:) = y; return
 116        end
 117        h = 2*h; lim = nstep; k = 1;
 118     else
 119        h = -h; lim = nstep+1; k = 0;
 120     end
 121
 122     for j=1:lim
 123         y = yi + h*pi;
 124         x(:) = y; fnew = dirn*feval(f,x,varargin{:}); nf = nf + 1;
 125         if fnew <= fmax, break, end
 126         fmax = fnew; k = k + 1;
 127         if fmax >= stopit(3)
 128            if trace
 129               fprintf('Comp. = %2.0f,  steps = %2.0f,  f = %9.4e*\n', i,j,fmax)
 130               fprintf('Exceeded target...quitting\n')
 131            end
 132            x(:) = y; return
 133         end
 134         h = 2*h;
 135    end
 136
 137    steps(i) = k;
 138    y = yi + 0.5*h*pi;
 139    if k == 0, y = yi; end
 140
 141    if trace
 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))
 144    end
 145
 146
 147    if nf >= stopit(2)
 148       if trace
 149          fprintf('Max no. of function evaluations exceeded...quitting\n')
 150       end
 151       x(:) = y; return
 152    end
 153
 154    if (fmax > flast && ~isempty (savit))
 155       x(:) = y;
 156       eval(['save ' savit ' x fmax nf'])
 157    end
 158
 159 end  % Loop over search directions.
 160
 161 if isequal(steps,zeros(n,1))
 162    if trace, fprintf('Stagnated...quitting\n'), end
 163    x(:) = y; return
 164 end
 165
 166 if fmax-fmax_old <= tol*abs(fmax_old)
 167    if trace, fprintf('Function values ''converged''...quitting\n'), end
 168    x(:) = y; return
 169 end
 170
 171 end %%%%%% Of outer loop.