Mercurial > matrix-functions

comparison matrixcomp/adsmax.m @ 0:8f23314345f4 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Create local repository for matrix toolboxes. Step #0 done.
author Antonio Pino Robles <data.script93@gmail.com>
date Wed, 06 May 2015 14:56:53 +0200
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:8f23314345f4
1 function [x, fmax, nf] = adsmax(f, x, stopit, savit, P, varargin)
2 %ADSMAX Alternating directions method for direct search optimization.
3 % [x, fmax, nf] = ADSMAX(FUN, x0, STOPIT, SAVIT, P) attempts to
4 % maximize the function FUN, using the starting vector x0.
5 % The alternating directions direct search method is used.
6 % Output arguments:
7 % x = vector yielding largest function value found,
8 % fmax = function value at x,
9 % nf = number of function evaluations.
10 % The iteration is terminated when either
11 % - the relative increase in function value between successive
12 % iterations is <= STOPIT(1) (default 1e-3),
13 % - STOPIT(2) function evaluations have been performed
14 % (default inf, i.e., no limit), or
15 % - a function value equals or exceeds STOPIT(3)
16 % (default inf, i.e., no test on function values).
17 % Progress of the iteration is not shown if STOPIT(5) = 0 (default 1).
18 % If a non-empty fourth parameter string SAVIT is present, then
19 % `SAVE SAVIT x fmax nf' is executed after each inner iteration.
20 % By default, the search directions are the co-ordinate directions.
21 % The columns of a fifth parameter matrix P specify alternative search
22 % directions (P = EYE is the default).
23 % NB: x0 can be a matrix. In the output argument, in SAVIT saves,
24 % and in function calls, x has the same shape as x0.
25 % ADSMAX(fun, x0, STOPIT, SAVIT, P, P1, P2,...) allows additional
26 % arguments to be passed to fun, via feval(fun,x,P1,P2,...).
27
28 % Reference:
29 % N. J. Higham, Optimization by direct search in matrix computations,
30 % SIAM J. Matrix Anal. Appl, 14(2): 317-333, 1993.
31 % N. J. Higham, Accuracy and Stability of Numerical Algorithms,
32 % Second edition, Society for Industrial and Applied Mathematics,
33 % Philadelphia, PA, 2002; sec. 20.5.
34
35 x0 = x(:); % Work with column vector internally.
36 n = length(x0);
37
38 mu = 1e-4; % Initial percentage change in components.
39 nstep = 25; % Max number of times to double or decrease h.
40
41 % Set up convergence parameters.
42 if nargin < 3 | isempty(stopit), stopit(1) = 1e-3; end
43 tol = stopit(1); % Required rel. increase in function value over one iteration.
44 if length(stopit) == 1, stopit(2) = inf; end % Max no. of f-evaluations.
45 if length(stopit) == 2, stopit(3) = inf; end % Default target for f-values.
46 if length(stopit) < 5, stopit(5) = 1; end % Default: show progress.
47 trace = stopit(5);
48 if nargin < 4, savit = []; end % File name for snapshots.
49
50 if nargin < 5 | isempty(P)
51 P = eye(n); % Matrix of search directions.
52 else
53 if ~isequal(size(P),[n n]) % Check for common error.
54 error('P must be of dimension the number of elements in x0.')
55 end
56 end
57
58 fmax = feval(f,x,varargin{:}); nf = 1;
59 if trace, fprintf('f(x0) = %9.4e\n', fmax), end
60
61 steps = zeros(n,1);
62 it = 0; y = x0;
63
64 while 1 % Outer loop.
65 it = it+1;
66 if trace, fprintf('Iter %2.0f (nf = %2.0f)\n', it, nf), end
67 fmax_old = fmax;
68
69 for i=1:n % Loop over search directions.
70
71 pi = P(:,i);
72 flast = fmax;
73 yi = y;
74 h = sign(pi'*yi)*norm(pi.*yi)*mu; % Initial step size.
75 if h == 0, h = max(norm(yi,inf),1)*mu; end
76 y = yi + h*pi;
77 x(:) = y; fnew = feval(f,x,varargin{:}); nf = nf + 1;
78 if fnew > fmax
79 fmax = fnew;
80 if fmax >= stopit(3)
81 if trace
82 fprintf('Comp. = %2.0f, steps = %2.0f, f = %9.4e*\n', i,0,fmax)
83 fprintf('Exceeded target...quitting\n')
84 end
85 x(:) = y; return
86 end
87 h = 2*h; lim = nstep; k = 1;
88 else
89 h = -h; lim = nstep+1; k = 0;
90 end
91
92 for j=1:lim
93 y = yi + h*pi;
94 x(:) = y; fnew = feval(f,x,varargin{:}); nf = nf + 1;
95 if fnew <= fmax, break, end
96 fmax = fnew; k = k + 1;
97 if fmax >= stopit(3)
98 if trace
99 fprintf('Comp. = %2.0f, steps = %2.0f, f = %9.4e*\n', i,j,fmax)
100 fprintf('Exceeded target...quitting\n')
101 end
102 x(:) = y; return
103 end
104 h = 2*h;
105 end
106
107 steps(i) = k;
108 y = yi + 0.5*h*pi;
109 if k == 0, y = yi; end
110
111 if trace
112 fprintf('Comp. = %2.0f, steps = %2.0f, f = %9.4e', i, k, fmax)
113 fprintf(' (%2.1f%%)\n', 100*(fmax-flast)/(abs(flast)+eps))
114 end
115
116
117 if nf >= stopit(2)
118 if trace
119 fprintf('Max no. of function evaluations exceeded...quitting\n')
120 end
121 x(:) = y; return
122 end
123
124 if fmax > flast & ~isempty(savit)
125 x(:) = y;
126 eval(['save ' savit ' x fmax nf'])
127 end
128
129 end % Loop over search directions.
130
131 if isequal(steps,zeros(n,1))
132 if trace, fprintf('Stagnated...quitting\n'), end
133 x(:) = y; return
134 end
135
136 if fmax-fmax_old <= tol*abs(fmax_old)
137 if trace, fprintf('Function values ''converged''...quitting\n'), end
138 x(:) = y; return
139 end
140
141 end %%%%%% Of outer loop.