simplex_optimization_routine

 Nelder-Mead like optimization routine.
 By default, we use standard values for the reflection, the expansion, the contraction and the shrink coefficients (alpha = 1, chi = 2, psi = 1 / 2 and σ = 1 / 2).
 See http://en.wikipedia.org/wiki/Nelder-Mead_method

 This routine uses the Nelder-Mead simplex (direct search) method.
 As chaining could reveal interesting to reach the solution neighborhood,
 the function automatically restarts from the current solution while
 amelioration is possible.

 INPUTS
 objective_function     [string]       Name of the iobjective function to be minimized.
 x                      [double]       n*1 vector, starting guess of the optimization routine.
 options                [structure]

 OUTPUTS

0001 function [x,fval,exitflag] = simplex_optimization_routine(objective_function,x,options,varargin)
0002 % Nelder-Mead like optimization routine.
0003 % By default, we use standard values for the reflection, the expansion, the contraction and the shrink coefficients (alpha = 1, chi = 2, psi = 1 / 2 and σ = 1 / 2).
0004 % See http://en.wikipedia.org/wiki/Nelder-Mead_method
0005 %
0006 % This routine uses the Nelder-Mead simplex (direct search) method.
0007 % As chaining could reveal interesting to reach the solution neighborhood,
0008 % the function automatically restarts from the current solution while
0009 % amelioration is possible.
0010 %
0011 % INPUTS
0012 % objective_function     [string]       Name of the iobjective function to be minimized.
0013 % x                      [double]       n*1 vector, starting guess of the optimization routine.
0014 % options                [structure]
0015 %
0016 % OUTPUTS
0017 %
0018 %
0019 
0020 % Copyright (C) 2010, 2011 Dynare Team
0021 %
0022 % This file is part of Dynare.
0023 %
0024 % Dynare is free software: you can redistribute it and/or modify
0025 % it under the terms of the GNU General Public License as published by
0026 % the Free Software Foundation, either version 3 of the License, or
0027 % (at your option) any later version.
0028 %
0029 % Dynare is distributed in the hope that it will be useful,
0030 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0031 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0032 % GNU General Public License for more details.
0033 %
0034 % You should have received a copy of the GNU General Public License
0035 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0036 global bayestopt_
0037 
0038 % Set verbose mode
0039 verbose = 2;
0040 
0041 % Set number of control variables.
0042 number_of_variables = length(x);
0043 
0044 % Set tolerance parameter.
0045 if isfield(options,'tolerance') && isfield(options.tolerance,'x')
0046     x_tolerance = options.tolerance.x;
0047 else
0048     x_tolerance = 1e-4;
0049 end
0050 
0051 % Set tolerance parameter.
0052 if isfield(options,'tolerance') && isfield(options.tolerance,'f')
0053     f_tolerance = options.tolerance.f;
0054 else
0055     f_tolerance = 1e-4;
0056 end
0057 
0058 % Set maximum number of iterations.
0059 if isfield(options,'maxiter')
0060     max_iterations = options.maxiter;
0061 else
0062     max_iterations = 1000;
0063 end
0064 % Set maximum number of iterations.
0065 if isfield(options,'maxfcall')
0066     max_func_calls = options.maxfcall;
0067 else
0068     max_func_calls = 500*number_of_variables;
0069 end
0070 
0071 % Set reflection parameter.
0072 if isfield(options,'reflection_parameter')
0073     if isfield(options.reflection_parameter,'value')
0074         rho = options.reflection_parameter.value;
0075     else
0076         rho = 1.0;
0077     end
0078     if isfield(options.reflection_parameter,'random')
0079         randomize_rho = options.reflection_parameter.random;
0080         lambda_rho = 1/rho;
0081     else
0082         randomize_rho = 0;
0083     end
0084 else
0085     rho = 1.0;
0086     randomize_rho = 0;
0087 end
0088 
0089 % Set expansion parameter.
0090 if isfield(options,'expansion_parameter')
0091     if isfield(options.expansion_parameter,'value')
0092         chi = options.expansion_parameter.value;
0093     else
0094         chi = 2.0;
0095     end
0096     if isfield(options.expansion_parameter,'random')
0097         randomize_chi = options.expansion_parameter.random;
0098         lambda_chi = 1/chi;
0099     else
0100         randomize_chi = 0;
0101     end
0102     if isfield(options.expansion_parameter,'optim')
0103         optimize_expansion_parameter = options.expansion_parameter.optim;
0104     else
0105         optimize_expansion_parameter = 0;
0106     end
0107 else
0108     chi = 2.0;
0109     randomize_chi = 0;
0110     optimize_expansion_parameter = 1;
0111 end
0112 
0113 % Set contraction parameter.
0114 if isfield(options,'contraction_parameter')
0115     if isfield(options.contraction_parameter,'value')
0116         psi = options.contraction_parameter.value;
0117     else
0118         psi = 0.5;
0119     end
0120     if isfield(options.contraction_parameter,'random')
0121         randomize_psi = options.expansion_parameter.random;
0122     else
0123         randomize_psi = 0;
0124     end
0125 else
0126     psi = 0.5;
0127     randomize_psi = 0;
0128 end
0129 
0130 % Set shrink parameter.
0131 if isfield(options,'shrink_parameter')
0132     if isfield(options.shrink_parameter,'value')
0133         sigma = options.shrink_parameter.value;
0134     else
0135         sigma = 0.5;
0136     end
0137     if isfield(options.shrink_parameter,'random')
0138         randomize_sigma = options.shrink_parameter.random;
0139     else
0140         randomize_sigma = 0;
0141     end
0142 else
0143     sigma = 0.5;
0144     randomize_sigma = 0;
0145 end
0146 
0147 % Set delta parameter.
0148 if isfield(options,'delta_parameter')
0149     delta = options.delta_parameter;
0150 else
0151     delta = 0.05;
0152 end
0153 DELTA = delta;
0154 zero_delta = delta/200;% To be used instead of delta if x(i) is zero.
0155 
0156 % Set max_no_improvements.
0157 if isfield(options,'max_no_improvements')
0158     max_no_improvements = options.max_no_improvements;
0159 else
0160     max_no_improvements = number_of_variables*10;
0161 end
0162 
0163 % Set vector of indices.
0164 unit_vector = ones(1,number_of_variables);
0165 trend_vector_1 = 1:number_of_variables;
0166 trend_vector_2 = 2:(number_of_variables+1);
0167 
0168 % Set initial simplex around the initial guess (x).
0169 if verbose
0170     for i=1:3, disp(' '), end
0171     disp('+----------------------+')
0172     disp(' SIMPLEX INITIALIZATION ')
0173     disp('+----------------------+')
0174     disp(' ')
0175 end
0176 initial_point = x;
0177 [initial_score,junk1,junk2,nopenalty] = feval(objective_function,x,varargin{:});
0178 if ~nopenalty
0179     error('simplex_optimization_routine:: Initial condition is wrong!')
0180 else
0181     [v,fv,delta] = simplex_initialization(objective_function,initial_point,initial_score,delta,1,varargin{:});
0182     func_count = number_of_variables + 1;
0183     iter_count = 1;
0184     if verbose
0185         disp(['Objective function value: ' num2str(fv(1))])
0186         disp(['Current parameter values: '])
0187         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');
0188         for i=1:number_of_variables
0189             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',bayestopt_.name{i},v(i,1), v(i,end), mean(v(i,:),2), min(v(i,:),[],2), max(v(i,:),[],2));
0190         end
0191         disp(' ')
0192     end
0193 end
0194 
0195 vold = v;
0196 
0197 no_improvements = 0;
0198 simplex_init = 1;
0199 simplex_iterations = 1;
0200 max_simplex_algo_iterations = 3;
0201 simplex_algo_iterations = 1;
0202 best_point = v(:,1);
0203 best_point_score = fv(1);
0204 convergence = 0;
0205 
0206 iter_no_improvement_break = 0;
0207 max_no_improvement_break = 1;
0208 
0209 while (func_count < max_func_calls) && (iter_count < max_iterations) && (simplex_algo_iterations<=max_simplex_algo_iterations)
0210     % Do we really need to continue?
0211     critF = max(abs(fv(1)-fv(trend_vector_2)));
0212     critX = max(max(abs(v(:,trend_vector_2)-v(:,unit_vector))));
0213     if critF <= max(f_tolerance,10*eps(fv(1))) && critX <= max(x_tolerance,10*eps(max(v(:,1))))
0214         convergence = 1;
0215     end
0216     % Set random reflection and expansion parameters if needed.
0217     if randomize_rho
0218         rho = -log(rand)/lambda_rho;
0219     end
0220     if randomize_chi
0221         chi = -log(rand)/lambda_chi;
0222     end
0223     % Set random contraction and shrink parameters if needed.
0224     if randomize_psi
0225         psi = rand;
0226     end
0227     if randomize_sigma
0228         sigma = rand;
0229     end
0230     % Compute the reflection point
0231     xbar = mean(v(:,trend_vector_1),2); % Average of the n best points.
0232     xr = xbar + rho*(xbar-v(:,end));
0233     x  = xr;
0234     fxr = feval(objective_function,x,varargin{:});
0235     func_count = func_count+1;
0236     if fxr < fv(1)% xr is better than previous best point v(:,1).
0237         % Calculate the expansion point
0238         xe = xbar + rho*chi*(xbar-v(:,end));
0239         x  = xe;
0240         fxe = feval(objective_function,x,varargin{:});
0241         func_count = func_count+1;
0242         if fxe < fxr% xe is even better than xr.
0243             if optimize_expansion_parameter
0244                 if verbose>1
0245                     disp('')
0246                     disp('')
0247                     disp('Compute optimal expansion...')
0248                 end
0249                 xee  = xbar + rho*chi*1.01*(xbar-v(:,end));
0250                 x    = xee;
0251                 fxee = feval(objective_function,x,varargin{:});
0252                 func_count = func_count+1;
0253                 if fxee<fxe
0254                     decrease = 1;
0255                     weight = rho*chi*1.02;
0256                     fxeee_old = fxee;
0257                     xeee_old  = xee;
0258                     if verbose>1
0259                         fprintf(1,'Weight =        ');
0260                     end
0261                     while decrease
0262                         weight = 1.02*weight;
0263                         if verbose>1
0264                             fprintf(1,'\b\b\b\b\b\b\b %6.4f',weight);
0265                         end
0266                         xeee   = xbar + weight*(xbar-v(:,end));
0267                         x      = xeee;
0268                         fxeee  = feval(objective_function,x,varargin{:});
0269                         func_count = func_count+1;
0270                         if (fxeee<fxeee_old) && -(fxeee-fxeee_old)>f_tolerance*10*fxeee_old
0271                             fxeee_old = fxeee;
0272                             xeee_old  = xeee;
0273                         else
0274                             decrease = 0;
0275                         end
0276                     end
0277                     if verbose>1
0278                         fprintf(1,'\n');
0279                     end
0280                     xe  = xeee_old;
0281                     fxe = fxeee_old;
0282                 else
0283                     decrease = 1;
0284                     weight = rho*chi;
0285                     fxeee_old = fxee;
0286                     xeee_old  = xee;
0287                     if verbose>1
0288                         fprintf(1,'Weight =        ');
0289                     end
0290                     while decrease
0291                         weight = weight/1.02;
0292                         if verbose>1
0293                             fprintf(1,'\b\b\b\b\b\b\b %6.4f',weight);
0294                         end
0295                         xeee   = xbar + weight*(xbar-v(:,end));
0296                         x      = xeee;
0297                         fxeee  = feval(objective_function,x,varargin{:});
0298                         func_count = func_count+1;
0299                         if (fxeee<fxeee_old) && -(fxeee-fxeee_old)>f_tolerance*10*fxeee_old
0300                             fxeee_old = fxeee;
0301                             xeee_old  = xeee;
0302                         else
0303                             decrease = 0;
0304                         end
0305                     end
0306                     if verbose>1
0307                         fprintf(1,'\n');
0308                     end
0309                     xe  = xeee_old;
0310                     fxe = fxeee_old;
0311                 end
0312                 if verbose>1
0313                     disp('Done!')
0314                     disp(' ')
0315                     disp(' ')
0316                 end
0317             end
0318             v(:,end) = xe;
0319             fv(end)  = fxe;
0320             move = 'expand';
0321         else% if xe is not better than xr.
0322             v(:,end) = xr;
0323             fv(end)  = fxr;
0324             move = 'reflect-1';
0325         end
0326     else% xr is not better than previous best point v(:,1).
0327         if fxr < fv(number_of_variables)% xr is better than previous point v(:,n).
0328             v(:,end) = xr;
0329             fv(end) = fxr;
0330             move = 'reflect-0';
0331         else% xr is not better than previous point v(:,n).
0332             if fxr < fv(end)% xr is better than previous worst point [=> outside contraction].
0333                 xc = (1 + psi*rho)*xbar - psi*rho*v(:,end);
0334                 x  = xc;
0335                 fxc = feval(objective_function,x,varargin{:});
0336                 func_count = func_count+1;
0337                 if fxc <= fxr
0338                     v(:,end) = xc;
0339                     fv(end) = fxc;
0340                     move = 'contract outside';
0341                 else
0342                     move = 'shrink';
0343                 end
0344             else% xr is the worst point [=> inside contraction].
0345                 xcc = (1-psi)*xbar + psi*v(:,end);
0346                 x   = xcc;
0347                 fxcc = feval(objective_function,x,varargin{:});
0348                 func_count = func_count+1;
0349                 if fxcc < fv(end)
0350                     v(:,end) = xcc;
0351                     fv(end) = fxcc;
0352                     move = 'contract inside';
0353                 else
0354                     % perform a shrink
0355                     move = 'shrink';
0356                 end
0357             end
0358             if strcmp(move,'shrink')
0359                 for j=trend_vector_2
0360                     v(:,j)=v(:,1)+sigma*(v(:,j) - v(:,1));
0361                     x = v(:,j);
0362                     fv(j) = feval(objective_function,x,varargin{:});
0363                 end
0364                 func_count = func_count + number_of_variables;
0365             end
0366         end
0367     end
0368     % Sort n+1 points by incresing order of the objective function values.
0369     [fv,sort_idx] = sort(fv);
0370     v = v(:,sort_idx);
0371     iter_count = iter_count + 1;
0372     simplex_iterations = simplex_iterations+1;
0373     if verbose>1
0374         disp(['Simplex iteration number: ' int2str(simplex_iterations) '-' int2str(simplex_init) '-' int2str(simplex_algo_iterations)])
0375         disp(['Simplex move:             ' move])
0376         disp(['Objective function value: ' num2str(fv(1))])
0377         disp(['Mode improvement:         ' num2str(best_point_score-fv(1))])
0378         disp(['Norm of dx:               ' num2str(norm(best_point-v(:,1)))])
0379         disp(['Norm of dSimplex:         ' num2str(norm(v-vold,'fro'))])
0380         disp(['Crit. f:                  ' num2str(critF)])
0381         disp(['Crit. x:                  ' num2str(critX)])
0382         disp(' ')
0383     end
0384     if verbose && max(abs(best_point-v(:,1)))>x_tolerance;
0385         if verbose<2
0386             disp(['Simplex iteration number: ' int2str(simplex_iterations) '-' int2str(simplex_init) '-' int2str(simplex_algo_iterations)])
0387             disp(['Objective function value: ' num2str(fv(1))])
0388             disp(['Mode improvement:         ' num2str(best_point_score-fv(1))])
0389             disp(['Norm of dx:               ' num2str(norm(best_point-v(:,1)))])
0390             disp(['Norm of dSimplex:         ' num2str(norm(v-vold,'fro'))])
0391             disp(['Crit. f:                  ' num2str(critF)])
0392             disp(['Crit. x:                  ' num2str(critX)])
0393             disp(' ')
0394         end
0395         disp(['Current parameter values: '])
0396         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');
0397         for i=1:number_of_variables
0398             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',bayestopt_.name{i}, v(i,1), v(i,end), mean(v(i,:),2), min(v(i,:),[],2), max(v(i,:),[],2));
0399         end
0400         disp(' ')
0401     end
0402     if abs(best_point_score-fv(1))<f_tolerance
0403         no_improvements = no_improvements+1;
0404     else
0405         no_improvements = 0;
0406     end
0407     best_point = v(:,1);
0408     best_point_score = fv(1);
0409     vold = v;
0410     if no_improvements>max_no_improvements
0411         if verbose
0412             disp(['NO SIGNIFICANT IMPROVEMENT AFTER ' int2str(no_improvements) ' ITERATIONS!'])
0413         end
0414         if simplex_algo_iterations<=max_simplex_algo_iterations
0415             % Compute the size of the simplex
0416             delta = delta*1.05;
0417             % Compute the new initial simplex.
0418             [v,fv,delta] = simplex_initialization(objective_function,best_point,best_point_score,delta,1,varargin{:});
0419             if verbose
0420                 disp(['(Re)Start with a lager simplex around the based on the best current '])
0421                 disp(['values for the control variables. '])
0422                 disp(['New value of delta (size of the new simplex) is: '])
0423                 for i=1:number_of_variables
0424                     fprintf(1,'%s: \t\t\t %+8.6f \n',bayestopt_.name{i}, delta(i));
0425                 end
0426             end
0427             % Reset counters
0428             no_improvements = 0;
0429             func_count = func_count + number_of_variables;
0430             iter_count = iter_count+1;
0431             iter_no_improvement_break = iter_no_improvement_break + 1;
0432             simplex_init = simplex_init+1;
0433             simplex_iterations = simplex_iterations+1;
0434             disp(' ')
0435             disp(' ')
0436         end
0437     end
0438     if ((func_count==max_func_calls) || (iter_count==max_iterations) || (iter_no_improvement_break==max_no_improvement_break) || convergence)
0439         [v,fv,delta] = simplex_initialization(objective_function,best_point,best_point_score,DELTA,1,varargin{:});
0440         if func_count==max_func_calls
0441             if verbose
0442                 disp(['MAXIMUM NUMBER OF OBJECTIVE FUNCTION CALLS EXCEEDED (' int2str(max_func_calls) ')!'])
0443             end
0444         elseif iter_count== max_iterations
0445             if verbose
0446                 disp(['MAXIMUM NUMBER OF ITERATIONS EXCEEDED (' int2str(max_iterations) ')!'])
0447             end
0448         elseif iter_no_improvement_break==max_no_improvement_break
0449             if verbose
0450                 disp(['MAXIMUM NUMBER OF SIMPLEX REINITIALIZATION EXCEEDED (' int2str(max_no_improvement_break) ')!'])
0451             end
0452             iter_no_improvement_break = 0;
0453             if simplex_algo_iterations==max_simplex_algo_iterations
0454                 max_no_improvements = Inf;% Do not stop until convergence is reached!
0455                 continue
0456             end
0457         else
0458             disp(['CONVERGENCE ACHIEVED AFTER ' int2str(simplex_iterations) ' ITERATIONS!'])
0459         end
0460         if simplex_algo_iterations<max_simplex_algo_iterations
0461             % Compute the size of the simplex
0462             delta = delta*1.05;
0463             % Compute the new initial simplex.
0464             [v,fv,delta] = simplex_initialization(objective_function,best_point,best_point_score,delta,1,varargin{:});
0465             if verbose
0466                 disp(['(Re)Start with a lager simplex around the based on the best current '])
0467                 disp(['values for the control variables. '])
0468                 disp(['New value of delta (size of the new simplex) is: '])
0469                 for i=1:number_of_variables
0470                     fprintf(1,'%s: \t\t\t %+8.6f \n',bayestopt_.name{i}, delta(i));
0471                 end
0472             end
0473             % Reset counters
0474             func_count=0;
0475             iter_count=0;
0476             convergence = 0;
0477             no_improvements = 0;
0478             func_count = func_count + number_of_variables;
0479             iter_count = iter_count+1;
0480             simplex_iterations = simplex_iterations+1;
0481             simplex_algo_iterations = simplex_algo_iterations+1;
0482             disp(' ')
0483             disp(' ')
0484         else
0485             break
0486         end
0487    end
0488 end% while loop.
0489 
0490 x(:) = v(:,1);
0491 fval = fv(:,1);
0492 exitflag = 1;
0493 
0494 if func_count>= max_func_calls
0495     disp('Maximum number of objective function calls has been exceeded!')
0496     exitflag = 0;
0497 end
0498 
0499 if iter_count>= max_iterations
0500     disp('Maximum number of iterations has been exceeded!')
0501     exitflag = 0;
0502 end
0503 
0504 
0505 
0506 
0507 function [v,fv,delta] = simplex_initialization(objective_function,point,point_score,delta,check_delta,varargin)
0508     n = length(point);
0509     v  = zeros(n,n+1);
0510     v(:,1) = point;
0511     fv = zeros(n+1,1);
0512     fv(1) = point_score;
0513     if length(delta)==1
0514         delta = repmat(delta,n,1);
0515     end
0516     for j = 1:n
0517         y = point;
0518         if y(j) ~= 0
0519             y(j) = (1 + delta(j))*y(j);
0520         else
0521             y(j) = zero_delta;
0522         end
0523         v(:,j+1) = y;
0524         x = y;
0525         [fv(j+1),junk1,junk2,nopenalty_flag] = feval(objective_function,x,varargin{:});
0526         if check_delta
0527             while ~nopenalty_flag
0528                 if y(j)~=0
0529                     delta(j) = delta(j)/1.1;
0530                 else
0531                     zero_delta = zero_delta/1.1;
0532                 end
0533                 y = point;
0534                 if y(j) ~= 0
0535                     y(j) = (1 + delta(j))*y(j);
0536                 else
0537                     y(j) = zero_delta;
0538                 end
0539                 v(:,j+1) = y;
0540                 x = y;
0541                 [fv(j+1),junk1,junk2,nopenalty_flag] = feval(objective_function,x,varargin{:});
0542             end
0543         end
0544     end
0545     % Sort by increasing order of the objective function values.
0546     [fv,sort_idx] = sort(fv);
0547     v = v(:,sort_idx);