simplex_optimization_routine.m 20.1 KB
 Johannes Pfeifer committed Jul 09, 2013 1 ``````function [x,fval,exitflag] = simplex_optimization_routine(objective_function,x,options,varargin) `````` Stéphane Adjemian committed Nov 16, 2013 2 3 `````` % Nelder-Mead like optimization routine (see http://en.wikipedia.org/wiki/Nelder-Mead_method) `````` Stéphane Adjemian committed May 10, 2011 4 ``````% `````` Stéphane Adjemian committed Nov 16, 2013 5 6 ``````% By default the standard values for the reflection, the expansion, the contraction % and the shrink coefficients are used (alpha = 1, chi = 2, psi = 1 / 2 and σ = 1 / 2). `````` Stéphane Adjemian committed May 10, 2011 7 ``````% `````` Stéphane Adjemian committed Nov 16, 2013 8 ``````% The routine automatically restarts from the current solution while amelioration is possible. `````` Stéphane Adjemian committed May 10, 2011 9 ``````% `````` Stéphane Adjemian committed Nov 16, 2013 10 11 12 13 14 15 ``````% INPUTS % o objective_function [string] Name of the objective function to be minimized. % o x [double] n*1 vector, starting guess of the optimization routine. % o options [structure] Options of this implementation of the simplex algorithm. % o varargin [cell of structures] Structures to be passed to the objective function: dataset_, % options_, M_, estim_params_, bayestopt_, and oo_. `````` Stéphane Adjemian committed May 10, 2011 16 ``````% `````` Stéphane Adjemian committed Nov 16, 2013 17 18 19 20 21 ``````% OUTPUTS % o x [double] n*1 vector, estimate of the optimal inputs. % o fval [double] scalar, value of the objective at the optimum. % o exitflag [integer] scalar equal to 0 or 1 (0 if the algorithm did not converge to % a minimum). `````` Stéphane Adjemian committed May 10, 2011 22 `````` `````` Johannes Pfeifer committed Jul 09, 2013 23 ``````% Copyright (C) 2010-2013 Dynare Team `````` Stéphane Adjemian committed May 10, 2011 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 ``````% % This file is part of Dynare. % % Dynare is free software: you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by % the Free Software Foundation, either version 3 of the License, or % (at your option) any later version. % % Dynare is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. % % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . % Set verbose mode verbose = 2; % Set number of control variables. number_of_variables = length(x); `````` Stéphane Adjemian committed Oct 08, 2013 46 ``````% get options. `````` Stéphane Adjemian committed Nov 16, 2013 47 48 ``````if isempty(options.maxfcall) max_func_calls = options.maxfcallfactor*number_of_variables `````` Stéphane Adjemian committed Oct 08, 2013 49 50 ``````end `````` Stéphane Adjemian committed May 10, 2011 51 ``````% Set tolerance parameter. `````` Stéphane Adjemian committed Mar 04, 2012 52 ``````if isfield(options,'tolerance') && isfield(options.tolerance,'x') `````` Stéphane Adjemian committed May 10, 2011 53 54 55 56 57 58 `````` x_tolerance = options.tolerance.x; else x_tolerance = 1e-4; end % Set tolerance parameter. `````` Stéphane Adjemian committed Mar 04, 2012 59 ``````if isfield(options,'tolerance') && isfield(options.tolerance,'f') `````` Stéphane Adjemian committed May 10, 2011 60 61 62 63 64 65 66 67 68 `````` f_tolerance = options.tolerance.f; else f_tolerance = 1e-4; end % Set maximum number of iterations. if isfield(options,'maxiter') max_iterations = options.maxiter; else `````` 69 `````` max_iterations = 5000; `````` Stéphane Adjemian committed May 10, 2011 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 ``````end % Set reflection parameter. if isfield(options,'reflection_parameter') if isfield(options.reflection_parameter,'value') rho = options.reflection_parameter.value; else rho = 1.0; end if isfield(options.reflection_parameter,'random') randomize_rho = options.reflection_parameter.random; lambda_rho = 1/rho; else randomize_rho = 0; end else rho = 1.0; randomize_rho = 0; end % Set expansion parameter. if isfield(options,'expansion_parameter') if isfield(options.expansion_parameter,'value') chi = options.expansion_parameter.value; else chi = 2.0; end if isfield(options.expansion_parameter,'random') randomize_chi = options.expansion_parameter.random; lambda_chi = 1/chi; else randomize_chi = 0; end if isfield(options.expansion_parameter,'optim') optimize_expansion_parameter = options.expansion_parameter.optim; else optimize_expansion_parameter = 0; end else chi = 2.0; randomize_chi = 0; optimize_expansion_parameter = 1; end % Set contraction parameter. if isfield(options,'contraction_parameter') if isfield(options.contraction_parameter,'value') psi = options.contraction_parameter.value; else psi = 0.5; end if isfield(options.contraction_parameter,'random') randomize_psi = options.expansion_parameter.random; else randomize_psi = 0; end else psi = 0.5; randomize_psi = 0; end % Set shrink parameter. if isfield(options,'shrink_parameter') if isfield(options.shrink_parameter,'value') sigma = options.shrink_parameter.value; else sigma = 0.5; end if isfield(options.shrink_parameter,'random') randomize_sigma = options.shrink_parameter.random; else randomize_sigma = 0; end else sigma = 0.5; randomize_sigma = 0; end % Set delta parameter. `````` Stéphane Adjemian committed Oct 08, 2013 149 ``````if isfield(options,'delta_parameter')% Size of the simplex `````` Stéphane Adjemian committed May 10, 2011 150 151 152 153 154 `````` delta = options.delta_parameter; else delta = 0.05; end DELTA = delta; `````` Marco Ratto committed Apr 29, 2012 155 ``````zero_delta = delta/200;% To be used instead of delta if x(i) is zero. `````` Stéphane Adjemian committed May 10, 2011 156 157 158 159 160 `````` % Set max_no_improvements. if isfield(options,'max_no_improvements') max_no_improvements = options.max_no_improvements; else `````` Stéphane Adjemian committed Mar 04, 2012 161 `````` max_no_improvements = number_of_variables*10; `````` Stéphane Adjemian committed May 10, 2011 162 163 164 165 166 167 168 169 170 ``````end % Set vector of indices. unit_vector = ones(1,number_of_variables); trend_vector_1 = 1:number_of_variables; trend_vector_2 = 2:(number_of_variables+1); % Set initial simplex around the initial guess (x). if verbose `````` Stéphane Adjemian committed Jul 10, 2013 171 `````` skipline(3) `````` Stéphane Adjemian committed May 10, 2011 172 173 174 `````` disp('+----------------------+') disp(' SIMPLEX INITIALIZATION ') disp('+----------------------+') `````` Stéphane Adjemian committed Jul 10, 2013 175 `````` skipline() `````` Stéphane Adjemian committed May 10, 2011 176 177 ``````end initial_point = x; `````` Marco Ratto committed Apr 29, 2012 178 ``````[initial_score,junk1,junk2,nopenalty] = feval(objective_function,x,varargin{:}); `````` Stéphane Adjemian committed May 10, 2011 179 180 181 ``````if ~nopenalty error('simplex_optimization_routine:: Initial condition is wrong!') else `````` Johannes Pfeifer committed Jul 09, 2013 182 `````` [v,fv,delta] = simplex_initialization(objective_function,initial_point,initial_score,delta,zero_delta,1,varargin{:}); `````` Stéphane Adjemian committed May 10, 2011 183 184 185 186 187 188 189 `````` func_count = number_of_variables + 1; iter_count = 1; if verbose disp(['Objective function value: ' num2str(fv(1))]) disp(['Current parameter values: ']) fprintf(1,'%s: \t\t\t %s \t\t\t %s \t\t\t %s \t\t\t %s \t\t\t %s \n','Names','Best point', 'Worst point', 'Mean values', 'Min values', 'Max values'); for i=1:number_of_variables `````` Stéphane Adjemian committed Nov 16, 2013 190 `````` fprintf(1,'%s: \t\t\t %+8.6f \t\t\t %+8.6f \t\t\t %+8.6f \t\t\t %+8.6f \t\t\t %+8.6f \n',varargin{5}.name{i},v(i,1), v(i,end), mean(v(i,:),2), min(v(i,:),[],2), max(v(i,:),[],2)); `````` Stéphane Adjemian committed May 10, 2011 191 `````` end `````` Stéphane Adjemian committed Jul 10, 2013 192 `````` skipline() `````` Stéphane Adjemian committed May 10, 2011 193 194 195 196 197 198 199 200 201 202 203 204 `````` end end vold = v; no_improvements = 0; simplex_init = 1; simplex_iterations = 1; max_simplex_algo_iterations = 3; simplex_algo_iterations = 1; best_point = v(:,1); best_point_score = fv(1); `````` Stéphane Adjemian committed Jun 24, 2013 205 `````` `````` Stéphane Adjemian committed May 10, 2011 206 ``````convergence = 0; `````` Stéphane Adjemian committed Jun 24, 2013 207 ``````tooslow = 0; `````` Stéphane Adjemian committed May 10, 2011 208 209 210 211 212 213 214 215 216 217 218 `````` iter_no_improvement_break = 0; max_no_improvement_break = 1; while (func_count < max_func_calls) && (iter_count < max_iterations) && (simplex_algo_iterations<=max_simplex_algo_iterations) % Do we really need to continue? critF = max(abs(fv(1)-fv(trend_vector_2))); critX = max(max(abs(v(:,trend_vector_2)-v(:,unit_vector)))); if critF <= max(f_tolerance,10*eps(fv(1))) && critX <= max(x_tolerance,10*eps(max(v(:,1)))) convergence = 1; end `````` Stéphane Adjemian committed Jun 24, 2013 219 220 221 `````` if critX <= x_tolerance^2 && critF>1 tooslow = 1; end `````` Stéphane Adjemian committed May 10, 2011 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 `````` % Set random reflection and expansion parameters if needed. if randomize_rho rho = -log(rand)/lambda_rho; end if randomize_chi chi = -log(rand)/lambda_chi; end % Set random contraction and shrink parameters if needed. if randomize_psi psi = rand; end if randomize_sigma sigma = rand; end % Compute the reflection point xbar = mean(v(:,trend_vector_1),2); % Average of the n best points. xr = xbar + rho*(xbar-v(:,end)); x = xr; fxr = feval(objective_function,x,varargin{:}); func_count = func_count+1; if fxr < fv(1)% xr is better than previous best point v(:,1). % Calculate the expansion point xe = xbar + rho*chi*(xbar-v(:,end)); x = xe; fxe = feval(objective_function,x,varargin{:}); func_count = func_count+1; if fxe < fxr% xe is even better than xr. if optimize_expansion_parameter if verbose>1 `````` Stéphane Adjemian committed Jul 10, 2013 251 `````` skipline(2) `````` Stéphane Adjemian committed May 10, 2011 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 `````` disp('Compute optimal expansion...') end xee = xbar + rho*chi*1.01*(xbar-v(:,end)); x = xee; fxee = feval(objective_function,x,varargin{:}); func_count = func_count+1; if fxee1 fprintf(1,'Weight = '); end while decrease weight = 1.02*weight; if verbose>1 fprintf(1,'\b\b\b\b\b\b\b %6.4f',weight); end xeee = xbar + weight*(xbar-v(:,end)); x = xeee; fxeee = feval(objective_function,x,varargin{:}); func_count = func_count+1; if (fxeeef_tolerance*10*fxeee_old fxeee_old = fxeee; xeee_old = xeee; else decrease = 0; end end if verbose>1 fprintf(1,'\n'); end xe = xeee_old; fxe = fxeee_old; else decrease = 1; weight = rho*chi; fxeee_old = fxee; xeee_old = xee; if verbose>1 fprintf(1,'Weight = '); end while decrease weight = weight/1.02; if verbose>1 fprintf(1,'\b\b\b\b\b\b\b %6.4f',weight); end xeee = xbar + weight*(xbar-v(:,end)); x = xeee; fxeee = feval(objective_function,x,varargin{:}); func_count = func_count+1; if (fxeeef_tolerance*10*fxeee_old fxeee_old = fxeee; xeee_old = xeee; else decrease = 0; end end if verbose>1 fprintf(1,'\n'); end xe = xeee_old; fxe = fxeee_old; end if verbose>1 disp('Done!') `````` Stéphane Adjemian committed Jul 10, 2013 319 `````` skipline(2) `````` Stéphane Adjemian committed May 10, 2011 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 `````` end end v(:,end) = xe; fv(end) = fxe; move = 'expand'; else% if xe is not better than xr. v(:,end) = xr; fv(end) = fxr; move = 'reflect-1'; end else% xr is not better than previous best point v(:,1). if fxr < fv(number_of_variables)% xr is better than previous point v(:,n). v(:,end) = xr; fv(end) = fxr; move = 'reflect-0'; else% xr is not better than previous point v(:,n). if fxr < fv(end)% xr is better than previous worst point [=> outside contraction]. xc = (1 + psi*rho)*xbar - psi*rho*v(:,end); `````` Stéphane Adjemian committed Mar 04, 2012 338 `````` x = xc; `````` Stéphane Adjemian committed May 10, 2011 339 340 341 342 343 344 345 346 347 348 349 `````` fxc = feval(objective_function,x,varargin{:}); func_count = func_count+1; if fxc <= fxr v(:,end) = xc; fv(end) = fxc; move = 'contract outside'; else move = 'shrink'; end else% xr is the worst point [=> inside contraction]. xcc = (1-psi)*xbar + psi*v(:,end); `````` Stéphane Adjemian committed Mar 04, 2012 350 `````` x = xcc; `````` Stéphane Adjemian committed May 10, 2011 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 `````` fxcc = feval(objective_function,x,varargin{:}); func_count = func_count+1; if fxcc < fv(end) v(:,end) = xcc; fv(end) = fxcc; move = 'contract inside'; else % perform a shrink move = 'shrink'; end end if strcmp(move,'shrink') for j=trend_vector_2 v(:,j)=v(:,1)+sigma*(v(:,j) - v(:,1)); x = v(:,j); fv(j) = feval(objective_function,x,varargin{:}); end func_count = func_count + number_of_variables; end end end % Sort n+1 points by incresing order of the objective function values. [fv,sort_idx] = sort(fv); v = v(:,sort_idx); iter_count = iter_count + 1; simplex_iterations = simplex_iterations+1; if verbose>1 disp(['Simplex iteration number: ' int2str(simplex_iterations) '-' int2str(simplex_init) '-' int2str(simplex_algo_iterations)]) disp(['Simplex move: ' move]) disp(['Objective function value: ' num2str(fv(1))]) disp(['Mode improvement: ' num2str(best_point_score-fv(1))]) disp(['Norm of dx: ' num2str(norm(best_point-v(:,1)))]) disp(['Norm of dSimplex: ' num2str(norm(v-vold,'fro'))]) disp(['Crit. f: ' num2str(critF)]) disp(['Crit. x: ' num2str(critX)]) `````` Stéphane Adjemian committed Jul 10, 2013 386 `````` skipline() `````` Stéphane Adjemian committed May 10, 2011 387 388 389 390 391 392 393 394 395 396 `````` end if verbose && max(abs(best_point-v(:,1)))>x_tolerance; if verbose<2 disp(['Simplex iteration number: ' int2str(simplex_iterations) '-' int2str(simplex_init) '-' int2str(simplex_algo_iterations)]) disp(['Objective function value: ' num2str(fv(1))]) disp(['Mode improvement: ' num2str(best_point_score-fv(1))]) disp(['Norm of dx: ' num2str(norm(best_point-v(:,1)))]) disp(['Norm of dSimplex: ' num2str(norm(v-vold,'fro'))]) disp(['Crit. f: ' num2str(critF)]) disp(['Crit. x: ' num2str(critX)]) `````` Stéphane Adjemian committed Jul 10, 2013 397 `````` skipline() `````` Stéphane Adjemian committed May 10, 2011 398 399 400 401 `````` end disp(['Current parameter values: ']) fprintf(1,'%s: \t\t\t %s \t\t\t %s \t\t\t %s \t\t\t %s \t\t\t %s \n','Names','Best point', 'Worst point', 'Mean values', 'Min values', 'Max values'); for i=1:number_of_variables `````` Stéphane Adjemian committed Nov 16, 2013 402 `````` fprintf(1,'%s: \t\t\t %+8.6f \t\t\t %+8.6f \t\t\t %+8.6f \t\t\t %+8.6f \t\t\t %+8.6f \n',varargin{5}.name{i}, v(i,1), v(i,end), mean(v(i,:),2), min(v(i,:),[],2), max(v(i,:),[],2)); `````` Stéphane Adjemian committed May 10, 2011 403 `````` end `````` Stéphane Adjemian committed Jul 10, 2013 404 `````` skipline() `````` Stéphane Adjemian committed May 10, 2011 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 `````` end if abs(best_point_score-fv(1))max_no_improvements if verbose disp(['NO SIGNIFICANT IMPROVEMENT AFTER ' int2str(no_improvements) ' ITERATIONS!']) end if simplex_algo_iterations<=max_simplex_algo_iterations % Compute the size of the simplex delta = delta*1.05; % Compute the new initial simplex. `````` Johannes Pfeifer committed Jul 09, 2013 422 `````` [v,fv,delta] = simplex_initialization(objective_function,best_point,best_point_score,delta,zero_delta,1,varargin{:}); `````` Stéphane Adjemian committed May 10, 2011 423 424 425 426 427 `````` if verbose disp(['(Re)Start with a lager simplex around the based on the best current ']) disp(['values for the control variables. ']) disp(['New value of delta (size of the new simplex) is: ']) for i=1:number_of_variables `````` Stéphane Adjemian committed Nov 16, 2013 428 `````` fprintf(1,'%s: \t\t\t %+8.6f \n',varargin{5}.name{i}, delta(i)); `````` Stéphane Adjemian committed May 10, 2011 429 430 431 432 433 434 435 436 437 `````` end end % Reset counters no_improvements = 0; func_count = func_count + number_of_variables; iter_count = iter_count+1; iter_no_improvement_break = iter_no_improvement_break + 1; simplex_init = simplex_init+1; simplex_iterations = simplex_iterations+1; `````` Stéphane Adjemian committed Jul 10, 2013 438 `````` skipline(2) `````` Stéphane Adjemian committed May 10, 2011 439 440 `````` end end `````` Stéphane Adjemian committed Jun 24, 2013 441 `````` if ((func_count==max_func_calls) || (iter_count==max_iterations) || (iter_no_improvement_break==max_no_improvement_break) || convergence || tooslow) `````` Johannes Pfeifer committed Jul 09, 2013 442 `````` [v,fv,delta] = simplex_initialization(objective_function,best_point,best_point_score,DELTA,zero_delta,1,varargin{:}); `````` Stéphane Adjemian committed May 10, 2011 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 `````` if func_count==max_func_calls if verbose disp(['MAXIMUM NUMBER OF OBJECTIVE FUNCTION CALLS EXCEEDED (' int2str(max_func_calls) ')!']) end elseif iter_count== max_iterations if verbose disp(['MAXIMUM NUMBER OF ITERATIONS EXCEEDED (' int2str(max_iterations) ')!']) end elseif iter_no_improvement_break==max_no_improvement_break if verbose disp(['MAXIMUM NUMBER OF SIMPLEX REINITIALIZATION EXCEEDED (' int2str(max_no_improvement_break) ')!']) end iter_no_improvement_break = 0; if simplex_algo_iterations==max_simplex_algo_iterations max_no_improvements = Inf;% Do not stop until convergence is reached! continue end `````` Stéphane Adjemian committed Jun 24, 2013 460 461 `````` elseif tooslow disp(['CONVERGENCE NOT ACHIEVED AFTER ' int2str(simplex_iterations) ' ITERATIONS! IMPROVING TOO SLOWLY!']) `````` Stéphane Adjemian committed May 10, 2011 462 463 464 465 466 467 468 `````` else disp(['CONVERGENCE ACHIEVED AFTER ' int2str(simplex_iterations) ' ITERATIONS!']) end if simplex_algo_iterations= max_func_calls disp('Maximum number of objective function calls has been exceeded!') exitflag = 0; end if iter_count>= max_iterations disp('Maximum number of iterations has been exceeded!') exitflag = 0; end `````` Johannes Pfeifer committed Jul 09, 2013 511 ``````function [v,fv,delta] = simplex_initialization(objective_function,point,point_score,delta,zero_delta,check_delta,varargin) `````` Stéphane Adjemian committed May 10, 2011 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 `````` n = length(point); v = zeros(n,n+1); v(:,1) = point; fv = zeros(n+1,1); fv(1) = point_score; if length(delta)==1 delta = repmat(delta,n,1); end for j = 1:n y = point; if y(j) ~= 0 y(j) = (1 + delta(j))*y(j); else y(j) = zero_delta; end v(:,j+1) = y; x = y; `````` Marco Ratto committed Apr 29, 2012 529 `````` [fv(j+1),junk1,junk2,nopenalty_flag] = feval(objective_function,x,varargin{:}); `````` Stéphane Adjemian committed May 10, 2011 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 `````` if check_delta while ~nopenalty_flag if y(j)~=0 delta(j) = delta(j)/1.1; else zero_delta = zero_delta/1.1; end y = point; if y(j) ~= 0 y(j) = (1 + delta(j))*y(j); else y(j) = zero_delta; end v(:,j+1) = y; x = y; `````` Marco Ratto committed Apr 29, 2012 545 `````` [fv(j+1),junk1,junk2,nopenalty_flag] = feval(objective_function,x,varargin{:}); `````` Stéphane Adjemian committed May 10, 2011 546 547 548 549 550 551 `````` end end end % Sort by increasing order of the objective function values. [fv,sort_idx] = sort(fv); v = v(:,sort_idx);``````