Mercurial > matrix-functions
comparison matrixcomp/nmsmax.m @ 0:8f23314345f4 draft
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] = nmsmax(fun, x, stopit, savit, varargin) | |
2 %NMSMAX Nelder-Mead simplex method for direct search optimization. | |
3 % [x, fmax, nf] = NMSMAX(FUN, x0, STOPIT, SAVIT) attempts to | |
4 % maximize the function FUN, using the starting vector x0. | |
5 % The Nelder-Mead 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 size of the simplex is <= STOPIT(1) | |
12 % (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 % The form of the initial simplex is determined by STOPIT(4): | |
18 % STOPIT(4) = 0: regular simplex (sides of equal length, the default) | |
19 % STOPIT(4) = 1: right-angled simplex. | |
20 % Progress of the iteration is not shown if STOPIT(5) = 0 (default 1). | |
21 % If a non-empty fourth parameter string SAVIT is present, then | |
22 % `SAVE SAVIT x fmax nf' is executed after each inner iteration. | |
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 % NMSMAX(fun, x0, STOPIT, SAVIT, P1, P2,...) allows additional | |
26 % arguments to be passed to fun, via feval(fun,x,P1,P2,...). | |
27 | |
28 % References: | |
29 % N. J. Higham, Optimization by direct search in matrix computations, | |
30 % SIAM J. Matrix Anal. Appl, 14(2): 317-333, 1993. | |
31 % C. T. Kelley, Iterative Methods for Optimization, Society for Industrial | |
32 % and Applied Mathematics, Philadelphia, PA, 1999. | |
33 | |
34 x0 = x(:); % Work with column vector internally. | |
35 n = length(x0); | |
36 | |
37 % Set up convergence parameters etc. | |
38 if nargin < 3 | isempty(stopit), stopit(1) = 1e-3; end | |
39 tol = stopit(1); % Tolerance for cgce test based on relative size of simplex. | |
40 if length(stopit) == 1, stopit(2) = inf; end % Max no. of f-evaluations. | |
41 if length(stopit) == 2, stopit(3) = inf; end % Default target for f-values. | |
42 if length(stopit) == 3, stopit(4) = 0; end % Default initial simplex. | |
43 if length(stopit) == 4, stopit(5) = 1; end % Default: show progress. | |
44 trace = stopit(5); | |
45 if nargin < 4, savit = []; end % File name for snapshots. | |
46 | |
47 V = [zeros(n,1) eye(n)]; | |
48 f = zeros(n+1,1); | |
49 V(:,1) = x0; f(1) = feval(fun,x,varargin{:}); | |
50 fmax_old = f(1); | |
51 | |
52 if trace, fprintf('f(x0) = %9.4e\n', f(1)), end | |
53 | |
54 k = 0; m = 0; | |
55 | |
56 % Set up initial simplex. | |
57 scale = max(norm(x0,inf),1); | |
58 if stopit(4) == 0 | |
59 % Regular simplex - all edges have same length. | |
60 % Generated from construction given in reference [18, pp. 80-81] of [1]. | |
61 alpha = scale / (n*sqrt(2)) * [ sqrt(n+1)-1+n sqrt(n+1)-1 ]; | |
62 V(:,2:n+1) = (x0 + alpha(2)*ones(n,1)) * ones(1,n); | |
63 for j=2:n+1 | |
64 V(j-1,j) = x0(j-1) + alpha(1); | |
65 x(:) = V(:,j); f(j) = feval(fun,x,varargin{:}); | |
66 end | |
67 else | |
68 % Right-angled simplex based on co-ordinate axes. | |
69 alpha = scale*ones(n+1,1); | |
70 for j=2:n+1 | |
71 V(:,j) = x0 + alpha(j)*V(:,j); | |
72 x(:) = V(:,j); f(j) = feval(fun,x,varargin{:}); | |
73 end | |
74 end | |
75 nf = n+1; | |
76 how = 'initial '; | |
77 | |
78 [temp,j] = sort(f); | |
79 j = j(n+1:-1:1); | |
80 f = f(j); V = V(:,j); | |
81 | |
82 alpha = 1; beta = 1/2; gamma = 2; | |
83 | |
84 while 1 %%%%%% Outer (and only) loop. | |
85 k = k+1; | |
86 | |
87 fmax = f(1); | |
88 if fmax > fmax_old | |
89 if ~isempty(savit) | |
90 x(:) = V(:,1); eval(['save ' savit ' x fmax nf']) | |
91 end | |
92 if trace | |
93 fprintf('Iter. %2.0f,', k) | |
94 fprintf([' how = ' how ' ']); | |
95 fprintf('nf = %3.0f, f = %9.4e (%2.1f%%)\n', nf, fmax, ... | |
96 100*(fmax-fmax_old)/(abs(fmax_old)+eps)) | |
97 end | |
98 end | |
99 fmax_old = fmax; | |
100 | |
101 %%% Three stopping tests from MDSMAX.M | |
102 | |
103 % Stopping Test 1 - f reached target value? | |
104 if fmax >= stopit(3) | |
105 msg = ['Exceeded target...quitting\n']; | |
106 break % Quit. | |
107 end | |
108 | |
109 % Stopping Test 2 - too many f-evals? | |
110 if nf >= stopit(2) | |
111 msg = ['Max no. of function evaluations exceeded...quitting\n']; | |
112 break % Quit. | |
113 end | |
114 | |
115 % Stopping Test 3 - converged? This is test (4.3) in [1]. | |
116 v1 = V(:,1); | |
117 size_simplex = norm(V(:,2:n+1)-v1(:,ones(1,n)),1) / max(1, norm(v1,1)); | |
118 if size_simplex <= tol | |
119 msg = sprintf('Simplex size %9.4e <= %9.4e...quitting\n', ... | |
120 size_simplex, tol); | |
121 break % Quit. | |
122 end | |
123 | |
124 % One step of the Nelder-Mead simplex algorithm | |
125 % NJH: Altered function calls and changed CNT to NF. | |
126 % Changed each `fr < f(1)' type test to `>' for maximization | |
127 % and re-ordered function values after sort. | |
128 | |
129 vbar = (sum(V(:,1:n)')/n)'; % Mean value | |
130 vr = (1 + alpha)*vbar - alpha*V(:,n+1); x(:) = vr; fr = feval(fun,x,varargin{:}); | |
131 nf = nf + 1; | |
132 vk = vr; fk = fr; how = 'reflect, '; | |
133 if fr > f(n) | |
134 if fr > f(1) | |
135 ve = gamma*vr + (1-gamma)*vbar; x(:) = ve; fe = feval(fun,x,varargin{:}); | |
136 nf = nf + 1; | |
137 if fe > f(1) | |
138 vk = ve; fk = fe; | |
139 how = 'expand, '; | |
140 end | |
141 end | |
142 else | |
143 vt = V(:,n+1); ft = f(n+1); | |
144 if fr > ft | |
145 vt = vr; ft = fr; | |
146 end | |
147 vc = beta*vt + (1-beta)*vbar; x(:) = vc; fc = feval(fun,x,varargin{:}); | |
148 nf = nf + 1; | |
149 if fc > f(n) | |
150 vk = vc; fk = fc; | |
151 how = 'contract,'; | |
152 else | |
153 for j = 2:n | |
154 V(:,j) = (V(:,1) + V(:,j))/2; | |
155 x(:) = V(:,j); f(j) = feval(fun,x,varargin{:}); | |
156 end | |
157 nf = nf + n-1; | |
158 vk = (V(:,1) + V(:,n+1))/2; x(:) = vk; fk = feval(fun,x,varargin{:}); | |
159 nf = nf + 1; | |
160 how = 'shrink, '; | |
161 end | |
162 end | |
163 V(:,n+1) = vk; | |
164 f(n+1) = fk; | |
165 [temp,j] = sort(f); | |
166 j = j(n+1:-1:1); | |
167 f = f(j); V = V(:,j); | |
168 | |
169 end %%%%%% End of outer (and only) loop. | |
170 | |
171 % Finished. | |
172 if trace, fprintf(msg), end | |
173 x(:) = V(:,1); |