stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

removed global objective_function_penalty_base 

added penalty_objective_function.m and penalty_hessian.m
time-shift
Michel Juillard 2015-09-20 15:28:26 +02:00
parent 9131472429
commit 5ade8d7c6f
13 changed files with 170 additions and 130 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

76
matlab/dsge_likelihood.m
View File

@ -134,8 +134,6 @@ function [fval,DLIK,Hess,exit_flag,SteadyState,trend_coeff,info,Model,DynareOpti
% AUTHOR(S) stephane DOT adjemian AT univ DASH lemans DOT FR
global objective_function_penalty_base
% Initialization of the returned variables and others...
fval = [];
SteadyState = [];
@ -181,9 +179,10 @@ end
% Return, with endogenous penalty, if some parameters are smaller than the lower bound of the prior domain.
if ~isequal(DynareOptions.mode_compute,1) && any(xparam1<BoundsInfo.lb)
k = find(xparam1<BoundsInfo.lb);
fval = objective_function_penalty_base+sum((BoundsInfo.lb(k)-xparam1(k)).^2);
fval = Inf;
exit_flag = 0;
info = 41;
info(1) = 41;
info(2) = sum((BoundsInfo.lb(k)-xparam1(k)).^2);
if analytic_derivation,
DLIK=ones(length(xparam1),1);
end
@ -193,9 +192,10 @@ end
% Return, with endogenous penalty, if some parameters are greater than the upper bound of the prior domain.
if ~isequal(DynareOptions.mode_compute,1) && any(xparam1>BoundsInfo.ub)
k = find(xparam1>BoundsInfo.ub);
fval = objective_function_penalty_base+sum((xparam1(k)-BoundsInfo.ub(k)).^2);
fval = Inf;
exit_flag = 0;
info = 42;
info(1) = 42;
info(2) = sum((xparam1(k)-BoundsInfo.ub(k)).^2);
if analytic_derivation,
DLIK=ones(length(xparam1),1);
end
@ -212,18 +212,19 @@ H = Model.H;
if ~issquare(Q) || EstimatedParameters.ncx || isfield(EstimatedParameters,'calibrated_covariances')
[Q_is_positive_definite, penalty] = ispd(Q);
if ~Q_is_positive_definite
fval = objective_function_penalty_base+penalty;
fval = Inf;
exit_flag = 0;
info = 43;
info(1) = 43;
info(2) = penalty;
return
end
if isfield(EstimatedParameters,'calibrated_covariances')
correct_flag=check_consistency_covariances(Q);
if ~correct_flag
penalty = sum(Q(EstimatedParameters.calibrated_covariances.position).^2);
fval = objective_function_penalty_base+penalty;
fval = Inf;
exit_flag = 0;
info = 71;
info(1) = 71;
info(2) = sum(Q(EstimatedParameters.calibrated_covariances.position).^2);
return
end
end
@ -233,18 +234,19 @@ end
if ~issquare(H) || EstimatedParameters.ncn || isfield(EstimatedParameters,'calibrated_covariances_ME')
[H_is_positive_definite, penalty] = ispd(H);
if ~H_is_positive_definite
fval = objective_function_penalty_base+penalty;
fval = Inf;
exit_flag = 0;
info = 44;
info(1) = 44;
info(2) = penalty;
return
end
if isfield(EstimatedParameters,'calibrated_covariances_ME')
correct_flag=check_consistency_covariances(H);
if ~correct_flag
penalty = sum(H(EstimatedParameters.calibrated_covariances_ME.position).^2);
fval = objective_function_penalty_base+penalty;
fval = Inf;
exit_flag = 0;
info = 72;
info(1) = 72;
info(2) = sum(H(EstimatedParameters.calibrated_covariances_ME.position).^2);
return
end
end
@ -261,16 +263,15 @@ end
% Return, with endogenous penalty when possible, if dynare_resolve issues an error code (defined in resol).
if info(1) == 1 || info(1) == 2 || info(1) == 5 || info(1) == 7 || info(1) == 8 || ...
info(1) == 22 || info(1) == 24 || info(1) == 19 || info(1) == 25 || info(1) == 10
fval = objective_function_penalty_base+1;
info = info(1);
fval = Inf;
info(2) = 1;
exit_flag = 0;
if analytic_derivation,
DLIK=ones(length(xparam1),1);
end
return
elseif info(1) == 3 || info(1) == 4 || info(1)==6 || info(1) == 20 || info(1) == 21 || info(1) == 23 || info(1)==26
fval = objective_function_penalty_base+info(2);
info = info(1);
fval = Inf;
exit_flag = 0;
if analytic_derivation,
DLIK=ones(length(xparam1),1);
@ -281,8 +282,7 @@ end
% check endogenous prior restrictions
info=endogenous_prior_restrictions(T,R,Model,DynareOptions,DynareResults);
if info(1),
fval = objective_function_penalty_base+info(2);
info = info(1);
fval = Inf;
exit_flag = 0;
if analytic_derivation,
DLIK=ones(length(xparam1),1);
@ -454,8 +454,9 @@ switch DynareOptions.lik_init
diffuse_periods = size(dlik,1);
end
if isnan(dLIK),
info = 45;
fval = objective_function_penalty_base + 100;
fval = Inf;
info(1) = 45;
info(2) = 100;
exit_flag = 0;
return
end
@ -683,8 +684,9 @@ if ((kalman_algo==1) || (kalman_algo==3))% Multivariate Kalman Filter
end
else
if isinf(LIK)
info = 66;
fval = objective_function_penalty_base+1;
fval = Inf;
info(1) = 66;
info(2) = 1;
exit_flag = 0;
return
end
@ -776,15 +778,17 @@ if analytic_derivation
end
if isnan(LIK)
info = 45;
fval = objective_function_penalty_base + 100;
fval = Inf;
info(1) = 45;
info(2) = 100;
exit_flag = 0;
return
end
if imag(LIK)~=0
info = 46;
fval = objective_function_penalty_base + 100;
fval = Inf;
info(1) = 46;
info(2) = 100;
exit_flag = 0;
return
end
@ -830,15 +834,17 @@ if DynareOptions.prior_restrictions.status
end
if isnan(fval)
info = 47;
fval = objective_function_penalty_base + 100;
fval = Inf;
info(1) = 47;
info(2) = 100;
exit_flag = 0;
return
end
if imag(fval)~=0
info = 48;
fval = objective_function_penalty_base + 100;
fval = Inf;
info(1) = 48;
info(2) = 100;
exit_flag = 0;
return
end
@ -846,7 +852,7 @@ end
% Update DynareOptions.kalman_algo.
DynareOptions.kalman_algo = kalman_algo;
if analytic_derivation==0 && nargout==2,
if analytic_derivation==0 && nargout > 1
lik=lik(start:end,:);
DLIK=[-lnprior; lik(:)];
end

43
matlab/dynare_estimation_1.m
View File

@ -238,26 +238,34 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
[xparam1, fval, exitflag, hh, options_, Scale] = dynare_minimize_objective(objective_function,xparam1,options_.mode_compute,options_,[bounds.lb bounds.ub],bayestopt_.name,bayestopt_,hh,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
fprintf('\nFinal value of minus the log posterior (or likelihood):%f \n', fval);
if isnumeric(options_.mode_compute) && options_.mode_compute==5 && options_.analytic_derivation==-1 %reset options changed by newrat
options_.analytic_derivation = options_analytic_derivation_old; %reset
elseif isnumeric(options_.mode_compute) && options_.mode_compute==6 %save scaling factor
if isequal(options_.mode_compute,5) && options_.analytic_derivation==-1 ...
% mode_compute = 5: reset options changed by newrat
options_.analytic_derivation = options_analytic_derivation_old; %reset
end
if isequal(options_.mode_compute,6) % mode_compute = 6 save scaling factor
save([M_.fname '_optimal_mh_scale_parameter.mat'],'Scale');
options_.mh_jscale = Scale;
bayestopt_.jscale = ones(length(xparam1),1)*Scale;
end
if ~isnumeric(options_.mode_compute) || ~isequal(options_.mode_compute,6) %always already computes covariance matrix
if options_.cova_compute == 1 %user did not request covariance not to be computed
if options_.analytic_derivation && strcmp(func2str(objective_function),'dsge_likelihood'),
ana_deriv_old = options_.analytic_derivation;
options_.analytic_derivation = 2;
[junk1, junk2, hh] = feval(objective_function,xparam1, ...
dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
options_.analytic_derivation = ana_deriv_old;
elseif ~isnumeric(options_.mode_compute) || ~(isequal(options_.mode_compute,5) && newratflag~=1),
if options_.cova_compute == 1 % user did not request covariance not to be computed
if options_.analytic_derivation && strcmp(func2str(objective_function),'dsge_likelihood')
options = options_.analytic_derivation;
options.analytic_derivation = 2;
[junk1, junk2, hh] = feval(objective_function,xparam1, ...
dataset_,dataset_info,options_,M_, ...
estim_params_,bayestopt_,bounds,oo_);
elseif isequal(options_.mode_compute,4) || ...
isequal(options_.mode_compute,5) % use penalty_hessian
% and fval as penalty base
if ~(isequal(options_.mode_compute,5) && newratflag ~= 1)
% with flag==0, we force to use the hessian from outer
% product gradient of optimizer 5
hh = reshape(hessian(objective_function,xparam1, ...
options_.gstep,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_),nx,nx);
hh = reshape(penalty_hessian(objective_function,xparam1, fval, ...
options_.gstep,dataset_, ...
dataset_info,options_,M_, ...
estim_params_,bayestopt_,bounds,oo_),nx,nx);
elseif isnumeric(options_.mode_compute) && isequal(options_.mode_compute,5)
% other numerical hessian options available with optimizer 5
%
@ -287,8 +295,15 @@ if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
end
options_.kalman_algo = kalman_algo0;
end
elseif isequal(options_.mode_compute,6) % mode_compute 6: hessian
% already computed, do
% nothing
else % all other cases: compute numerical hessian
hh = reshape(hessian(objective_function,xparam1, ...
options_.gstep,dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_),nx,nx);
end
end
parameter_names = bayestopt_.name;
if options_.cova_compute || options_.mode_compute==5 || options_.mode_compute==6
save([M_.fname '_mode.mat'],'xparam1','hh','parameter_names','fval');

5
matlab/optimization/csminit1.m
View File

@ -1,9 +1,10 @@
function [fhat,xhat,fcount,retcode] = csminit1(fcn,x0,f0,g0,badg,H0,Verbose,varargin)
function [fhat,xhat,fcount,retcode] = csminit1(fcn,x0,penalty,f0,g0,badg,H0,Verbose,varargin)
% [fhat,xhat,fcount,retcode] = csminit1(fcn,x0,f0,g0,badg,H0,varargin)
%
% Inputs:
% fcn: [string] string naming the objective function to be minimized
% x0: [npar by 1] initial value of the parameter vector
% penalty: [sca;ar] variable penalty in case of failure of objective function
% g0: [npar by 1] initial value of the gradient vector
% H0: [npar by npar] initial value for the inverse Hessian. Must be positive definite.
% varargin: Optional additional inputs that get handed off to fcn each
@ -137,7 +138,7 @@ else
dxtest=x0+dx*lambda;
end
% home
f = feval(fcn,dxtest,varargin{:});
f = penalty_objective_function(dxtest,fcn,penalty,varargin{:});
%ARGLIST
%f = feval(fcn,dxtest,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13);
% f = feval(fcn,x0+dx*lambda,P1,P2,P3,P4,P5,P6,P7,P8);

62
matlab/optimization/csminwel1.m
View File

@ -59,7 +59,8 @@ function [fh,xh,gh,H,itct,fcount,retcodeh] = csminwel1(fcn,x0,H0,grad,crit,nit,m
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global objective_function_penalty_base
% initialize variable penalty
penalty = 1e8;
fh = [];
xh = [];
@ -72,6 +73,15 @@ fcount=0;
gh = [];
H = [];
retcodeh = [];
% force fcn, grad to function handle
if ischar(fcn)
fcn = str2func(fcn);
end
if ischar(grad)
grad = str2func(grad)
end
%tailstr = ')';
%stailstr = [];
% Lines below make the number of Pi's optional. This is inefficient, though, and precludes
@ -83,7 +93,7 @@ retcodeh = [];
% stailstr=[' P' num2str(i) stailstr];
%end
[f0,junk1,junk2,cost_flag] = feval(fcn,x0,varargin{:});
[f0,junk1,junk2,cost_flag] = penalty_objective_function(x0,fcn,penalty,varargin{:});
if ~cost_flag
disp_verbose('Bad initial parameter.',Verbose)
@ -91,9 +101,9 @@ if ~cost_flag
end
if NumGrad
[g, badg]=get_num_grad(method,fcn,f0,x0,epsilon,varargin{:});
[g, badg]=get_num_grad(method,fcn,penalty,f0,x0,epsilon,varargin{:});
elseif ischar(grad)
[g,badg] = feval(grad,x0,varargin{:});
[g,badg] = grad(x0,varargin{:});
else
g=junk1;
badg=0;
@ -105,7 +115,7 @@ H=H0;
cliff=0;
while ~done
% penalty for dsge_likelihood and dsge_var_likelihood
objective_function_penalty_base = f;
penalty = f;
g1=[]; g2=[]; g3=[];
%addition fj. 7/6/94 for control
@ -115,7 +125,7 @@ while ~done
% disp_verbose([sprintf('x = ') sprintf('%15.8g %15.8g %15.8g %15.8g\n',x)]);
%-------------------------
itct=itct+1;
[f1, x1, fc, retcode1] = csminit1(fcn,x,f,g,badg,H,Verbose,varargin{:});
[f1, x1, fc, retcode1] = csminit1(fcn,x,penalty,f,g,badg,H,Verbose,varargin{:});
fcount = fcount+fc;
% erased on 8/4/94
% if (retcode == 1) || (abs(f1-f) < crit)
@ -130,11 +140,11 @@ while ~done
wall1=1; badg1=1;
else
if NumGrad
[g1, badg1]=get_num_grad(method,fcn,f1,x1,epsilon,varargin{:});
[g1, badg1]=get_num_grad(method,fcn,penalty,f1,x1,epsilon,varargin{:});
elseif ischar(grad),
[g1, badg1] = feval(grad,x1,varargin{:});
[g1, badg1] = grad(x1,varargin{:});
else
[junk1,g1,junk2, cost_flag] = feval(fcn,x1,varargin{:});
[junk1,g1,junk2, cost_flag] = penalty_objective_function(x1,fcn,penalty,varargin{:});
badg1 = ~cost_flag;
end
wall1=badg1;
@ -150,18 +160,18 @@ while ~done
%fcliff=fh;xcliff=xh;
Hcliff=H+diag(diag(H).*rand(nx,1));
disp_verbose('Cliff. Perturbing search direction.',Verbose)
[f2, x2, fc, retcode2] = csminit1(fcn,x,f,g,badg,Hcliff,Verbose,varargin{:});
[f2, x2, fc, retcode2] = csminit1(fcn,x,penalty,f,g,badg,Hcliff,Verbose,varargin{:});
fcount = fcount+fc; % put by Jinill
if f2 < f
if retcode2==2 || retcode2==4
wall2=1; badg2=1;
else
if NumGrad
[g2, badg2]=get_num_grad(method,fcn,f2,x2,epsilon,varargin{:});
[g2, badg2]=get_num_grad(method,fcn,penalty,f2,x2,epsilon,varargin{:});
elseif ischar(grad),
[g2, badg2] = feval(grad,x2,varargin{:});
[g2, badg2] = grad(x2,varargin{:});
else
[junk1,g2,junk2, cost_flag] = feval(fcn,x1,varargin{:});
[junk1,g2,junk2, cost_flag] = penalty_objective_function(x1,fcn,penalty,varargin{:});
badg2 = ~cost_flag;
end
wall2=badg2;
@ -182,18 +192,18 @@ while ~done
if(size(x0,2)>1)
gcliff=gcliff';
end
[f3, x3, fc, retcode3] = csminit1(fcn,x,f,gcliff,0,eye(nx),Verbose,varargin{:});
[f3, x3, fc, retcode3] = csminit1(fcn,x,penalty,f,gcliff,0,eye(nx),Verbose,varargin{:});
fcount = fcount+fc; % put by Jinill
if retcode3==2 || retcode3==4
wall3=1;
badg3=1;
else
if NumGrad
[g3, badg3]=get_num_grad(method,fcn,f3,x3,epsilon,varargin{:});
[g3, badg3]=get_num_grad(method,fcn,penalty,f3,x3,epsilon,varargin{:});
elseif ischar(grad),
[g3, badg3] = feval(grad,x3,varargin{:});
[g3, badg3] = grad(x3,varargin{:});
else
[junk1,g3,junk2, cost_flag] = feval(fcn,x1,varargin{:});
[junk1,g3,junk2, cost_flag] = penalty_objective_function(x1,fcn,penalty,varargin{:});
badg3 = ~cost_flag;
end
wall3=badg3;
@ -249,11 +259,11 @@ while ~done
end
if nogh
if NumGrad
[gh, badgh]=get_num_grad(method,fcn,fh,xh,epsilon,varargin{:});
[gh, badgh]=get_num_grad(method,fcn,penalty,fh,xh,epsilon,varargin{:});
elseif ischar(grad),
[gh, badgh] = feval(grad, xh,varargin{:});
[gh, badgh] = grad(xh,varargin{:});
else
[junk1,gh,junk2, cost_flag] = feval(fcn,x1,varargin{:});
[junk1,gh,junk2, cost_flag] = penalty_objective_function(x1,penalty,varargin{:});
badgh = ~cost_flag;
end
end
@ -303,18 +313,18 @@ end
end
function [g, badg]=get_num_grad(method,fcn,f0,x0,epsilon,varargin)
function [g, badg]=get_num_grad(method,fcn,penalty,f0,x0,epsilon,varargin)
switch method
case 2
[g,badg] = numgrad2(fcn, f0, x0, epsilon, varargin{:});
[g,badg] = numgrad2(fcn, f0, x0, penalty, epsilon, varargin{:});
case 3
[g,badg] = numgrad3(fcn, f0, x0, epsilon, varargin{:});
[g,badg] = numgrad3(fcn, f0, x0, penalty, epsilon, varargin{:});
case 5
[g,badg] = numgrad5(fcn, f0, x0, epsilon, varargin{:});
[g,badg] = numgrad5(fcn, f0, x0, penalty, epsilon, varargin{:});
case 13
[g,badg] = numgrad3_(fcn, f0, x0, epsilon, varargin{:});
[g,badg] = numgrad3_(fcn, f0, x0, penalty, epsilon, varargin{:});
case 15
[g,badg] = numgrad5_(fcn, f0, x0, epsilon, varargin{:});
[g,badg] = numgrad5_(fcn, f0, x0, penalty, epsilon, varargin{:});
otherwise
error('csminwel1: Unknown method for gradient evaluation!')
end

3
matlab/optimization/dynare_minimize_objective.m
View File

@ -244,7 +244,7 @@ switch minimizer_algorithm
Save_files = 0;
Verbose = 0;
end
[opt_par_values,hessian_mat,gg,fval,invhess] = newrat(objective_function,start_par_value,analytic_grad,crit,nit,0,Verbose, Save_files,varargin{:});
[opt_par_values,hessian_mat,gg,fval,invhess] = newrat(objective_function,start_par_value,analytic_grad,crit,nit,newratflag,Verbose, Save_files,varargin{:});
%hessian_mat is the plain outer product gradient Hessian
case 6
[opt_par_values, hessian_mat, Scale, fval] = gmhmaxlik(objective_function, start_par_value, ...
@ -450,3 +450,4 @@ function [LB, UB]=set_bounds_to_finite_values(bounds, huge_number)
UB=bounds(:,2);
UB(isinf(UB))=huge_number;
end

10
matlab/optimization/mr_gstep.m
View File

@ -1,4 +1,4 @@
function [f0, x, ig] = mr_gstep(h1,x,func0,htol0,Verbose,Save_files,varargin)
function [f0, x, ig] = mr_gstep(h1,x,func0,penalty,htol0,Verbose,Save_files,varargin)
% function [f0, x, ig] = mr_gstep(h1,x,func0,htol0,varargin)
%
% Gibbs type step in optimisation
@ -39,7 +39,7 @@ if isempty(htol0)
else
htol = htol0;
end
f0=feval(func0,x,varargin{:});
f0=penalty_objective_function(x,func0,penalty,varargin{:});
xh1=x;
f1=zeros(size(f0,1),n);
@ -53,10 +53,10 @@ while i<n
hcheck=0;
dx=[];
xh1(i)=x(i)+h1(i);
fx = feval(func0,xh1,varargin{:});
fx = penalty_objective_function(xh1,func0,penalty,varargin{:});
f1(:,i)=fx;
xh1(i)=x(i)-h1(i);
fx = feval(func0,xh1,varargin{:});
fx = penalty_objective_function(xh1,func0,penalty,varargin{:});
f_1(:,i)=fx;
if hcheck && htol<1
htol=min(1,max(min(abs(dx))*2,htol*10));
@ -69,7 +69,7 @@ while i<n
gg(i)=(f1(i)'-f_1(i)')./(2.*h1(i));
hh(i) = 1/max(1.e-9,abs( (f1(i)+f_1(i)-2*f0)./(h1(i)*h1(i)) ));
if gg(i)*(hh(i)*gg(i))/2 > htol
[f0 x fc retcode] = csminit1(func0,x,f0,gg,0,diag(hh),Verbose,varargin{:});
[f0 x fc retcode] = csminit1(func0,x,penalty,f0,gg,0,diag(hh),Verbose,varargin{:});
ig(i)=1;
if Verbose
fprintf(['Done for param %s = %8.4f\n'],varargin{6}.name{i},x(i))

21
matlab/optimization/mr_hessian.m
View File

@ -1,5 +1,5 @@
function [hessian_mat, gg, htol1, ihh, hh_mat0, hh1] = mr_hessian(init,x,func,hflag,htol0,varargin)
% [hessian_mat, gg, htol1, ihh, hh_mat0, hh1] = mr_hessian(init,x,func,hflag,htol0,varargin)
function [hessian_mat, gg, htol1, ihh, hh_mat0, hh1] = mr_hessian(init,x,func,penalty,hflag,htol0,varargin)
% [hessian_mat, gg, htol1, ihh, hh_mat0, hh1] = mr_hessian(init,x,func,penaltyhflag,htol0,varargin)
%
% numerical gradient and Hessian, with 'automatic' check of numerical
% error
@ -10,6 +10,7 @@ function [hessian_mat, gg, htol1, ihh, hh_mat0, hh1] = mr_hessian(init,x,func,hf
% - the log-likelihood AND the single contributions at times t=1,...,T
% of the log-likelihood to compute outer product gradient
% x = parameter values
% penalty = endogenous penalty in case of objective function failure
% hflag = 0, Hessian computed with outer product gradient, one point
% increments for partial derivatives in gradients
% hflag = 1, 'mixed' Hessian: diagonal elements computed with numerical second order derivatives
@ -60,7 +61,7 @@ if init
return
end
[f0, ff0]=feval(func,x,varargin{:});
[f0, ff0]=penalty_objective_function(x,func,penalty,varargin{:});
h2=varargin{7}.ub-varargin{7}.lb;
hmax=varargin{7}.ub-x;
hmax=min(hmax,x-varargin{7}.lb);
@ -92,7 +93,7 @@ while i<n
hcheck=0;
xh1(i)=x(i)+h1(i);
try
[fx, ffx]=feval(func,xh1,varargin{:});
[fx, ffx]=penalty_objective_function(xh1,func,penalty,varargin{:});
catch
fx=1.e8;
end
@ -113,7 +114,7 @@ while i<n
h1(i) = max(h1(i),1.e-10);
xh1(i)=x(i)+h1(i);
try
[fx, ffx]=feval(func,xh1,varargin{:});
[fx, ffx]=penalty_objective_function(xh1,func,penalty,varargin{:});
catch
fx=1.e8;
end
@ -122,14 +123,14 @@ while i<n
h1(i)= htol/abs(dx(it))*h1(i);
xh1(i)=x(i)+h1(i);
try
[fx, ffx]=feval(func,xh1,varargin{:});
[fx, ffx]=penalty_objective_function(xh1,func,penalty,varargin{:});
catch
fx=1.e8;
end
while (fx-f0)==0
h1(i)= h1(i)*2;
xh1(i)=x(i)+h1(i);
[fx, ffx]=feval(func,xh1,varargin{:});
[fx, ffx]=penalty_objective_function(xh1,func,penalty,varargin{:});
ic=1;
end
end
@ -150,7 +151,7 @@ while i<n
end
end
xh1(i)=x(i)-h1(i);
[fx, ffx]=feval(func,xh1,varargin{:});
[fx, ffx]=penalty_objective_function(xh1,func,penalty,varargin{:});
f_1(:,i)=fx;
if outer_product_gradient,
if any(isnan(ffx)) || isempty(ffx),
@ -190,8 +191,8 @@ if outer_product_gradient,
xh1(j)=x(j)+h_1(j);
xh_1(i)=x(i)-h1(i);
xh_1(j)=x(j)-h_1(j);
temp1 = feval(func,xh1,varargin{:});
temp2 = feval(func,xh_1,varargin{:});
temp1 = penalty_objective_function(xh1,func,penalty,varargin{:});
temp2 = penalty_objective_function(xh,func,penalty_1,varargin{:});
hessian_mat(:,(i-1)*n+j)=-(-temp1 -temp2+temp(:,i)+temp(:,j))./(2*h1(i)*h_1(j));
xh1(i)=x(i);
xh1(j)=x(j);

44
matlab/optimization/newrat.m
View File

@ -49,6 +49,7 @@ function [xparam1, hh, gg, fval, igg] = newrat(func0, x, analytic_derivation, ft
global objective_function_penalty_base
penalty = 1e8;
icount=0;
nx=length(x);
@ -62,17 +63,22 @@ htol=htol_base;
htol0=htol_base;
gibbstol=length(varargin{6}.pshape)/50; %25;
% force fcn, grad to function handle
if ischar(func0)
func0 = str2func(func0);
end
% func0 = str2func([func2str(func0),'_hh']);
% func0 = func0;
[fval0,gg,hh]=feval(func0,x,varargin{:});
[fval0,gg,hh]=penalty_objective_function(x,func0,penalty,varargin{:});
fval=fval0;
% initialize mr_gstep and mr_hessian
outer_product_gradient=1;
if isempty(hh)
mr_hessian(1,x,[],[],[],varargin{:});
[dum, gg, htol0, igg, hhg, h1]=mr_hessian(0,x,func0,flagit,htol,varargin{:});
mr_hessian(1,x,[],[],[],[],varargin{:});
[dum, gg, htol0, igg, hhg, h1]=mr_hessian(0,x,func0,penalty,flagit,htol,varargin{:});
if isempty(dum),
outer_product_gradient=0;
igg = 1e-4*eye(nx);
@ -121,12 +127,12 @@ while norm(gg)>gtol && check==0 && jit<nit
jit=jit+1;
tic
icount=icount+1;
objective_function_penalty_base = fval0(icount);
penalty = fval0(icount);
disp_verbose([' '],Verbose)
disp_verbose(['Iteration ',num2str(icount)],Verbose)
[fval,x0,fc,retcode] = csminit1(func0,xparam1,fval0(icount),gg,0,H,Verbose,varargin{:});
[fval,x0,fc,retcode] = csminit1(func0,xparam1,penalty,fval0(icount),gg,0,H,Verbose,varargin{:});
if igrad
[fval1,x01,fc,retcode1] = csminit1(func0,x0,fval,gg,0,inx,Verbose,varargin{:});
[fval1,x01,fc,retcode1] = csminit1(func0,x0,penalty,fval,gg,0,inx,Verbose,varargin{:});
if (fval-fval1)>1
disp_verbose('Gradient step!!',Verbose)
else
@ -138,23 +144,23 @@ while norm(gg)>gtol && check==0 && jit<nit
if length(find(ig))<nx
ggx=ggx*0;
ggx(find(ig))=gg(find(ig));
if analytic_derivation,
hhx=hh;
else
hhx = reshape(dum,nx,nx);
end
% if analytic_derivation,
hhx=hh;
% else
% hhx = reshape(dum,nx,nx);
% end
iggx=eye(length(gg));
iggx(find(ig),find(ig)) = inv( hhx(find(ig),find(ig)) );
[fvala,x0,fc,retcode] = csminit1(func0,x0,fval,ggx,0,iggx,Verbose,varargin{:});
[fvala,x0,fc,retcode] = csminit1(func0,x0,penalty,fval,ggx,0,iggx,Verbose,varargin{:});
end
[fvala, x0, ig] = mr_gstep(h1,x0,func0,htol,Verbose,Save_files,varargin{:});
[fvala, x0, ig] = mr_gstep(h1,x0,func0,penalty,htol,Verbose,Save_files,varargin{:});
nig=[nig ig];
disp_verbose('Sequence of univariate steps!!',Verbose)
fval=fvala;
if (fval0(icount)-fval)<ftol && flagit==0
disp_verbose('Try diagonal Hessian',Verbose)
ihh=diag(1./(diag(hhg)));
[fval2,x0,fc,retcode2] = csminit1(func0,x0,fval,gg,0,ihh,Verbose,varargin{:});
[fval2,x0,fc,retcode2] = csminit1(func0,x0,penalty,fval,gg,0,ihh,Verbose,varargin{:});
if (fval-fval2)>=ftol
disp_verbose('Diagonal Hessian successful',Verbose)
end
@ -163,7 +169,7 @@ while norm(gg)>gtol && check==0 && jit<nit
if (fval0(icount)-fval)<ftol && flagit==0
disp_verbose('Try gradient direction',Verbose)
ihh0=inx.*1.e-4;
[fval3,x0,fc,retcode3] = csminit1(func0,x0,fval,gg,0,ihh0,Verbose,varargin{:});
[fval3,x0,fc,retcode3] = csminit1(func0,x0,penalty,fval,gg,0,ihh0,Verbose,varargin{:});
if (fval-fval3)>=ftol
disp_verbose('Gradient direction successful',Verbose)
end
@ -176,14 +182,14 @@ while norm(gg)>gtol && check==0 && jit<nit
disp_verbose('No further improvement is possible!',Verbose)
check=1;
if analytic_derivation,
[fvalx,gg,hh]=feval(func0,xparam1,varargin{:});
[fvalx,gg,hh]=penalty_objective_function(xparam1,func0,varargin{:});
hhg=hh;
H = inv(hh);
else
if flagit==2
hh=hh0;
elseif flagg>0
[dum, gg, htol0, igg, hhg,h1]=mr_hessian(0,xparam1,func0,flagg,ftol0,varargin{:});
[dum, gg, htol0, igg, hhg,h1]=mr_hessian(0,xparam1,func0,penalty,flagg,ftol0,varargin{:});
if flagg==2
hh = reshape(dum,nx,nx);
ee=eig(hh);
@ -223,7 +229,7 @@ while norm(gg)>gtol && check==0 && jit<nit
save m1.mat x fval0 nig
end
end
[dum, gg, htol0, igg, hhg, h1]=mr_hessian(0,xparam1,func0,flagit,htol,varargin{:});
[dum, gg, htol0, igg, hhg, h1]=mr_hessian(0,xparam1,func0,penalty,flagit,htol,varargin{:});
if isempty(dum),
outer_product_gradient=0;
end
@ -252,7 +258,7 @@ while norm(gg)>gtol && check==0 && jit<nit
H = igg;
end
elseif analytic_derivation,
[fvalx,gg,hh]=feval(func0,xparam1,varargin{:});
[fvalx,gg,hh]=penalty_objective_function(xparam1,func0,varargin{:});
hhg=hh;
H = inv(hh);
end

6
matlab/optimization/numgrad2.m
View File

@ -1,4 +1,4 @@
function [g, badg] = numgrad2(fcn,f0,x,epsilon,varargin)
function [g, badg] = numgrad2(fcn,f0,x,penalty,epsilon,varargin)
% Original file downloaded from:
% http://sims.princeton.edu/yftp/optimize/mfiles/numgrad.m
@ -33,11 +33,11 @@ for i=1:n
xiold = x(i);
if right_derivative
x(i) = x(i)+h;
fh = feval(fcn, x, varargin{:});
fh = penalty_objective_function(x, fcn, penalty, varargin{:});
g0 = (fh-f0)/h;
else
x(i) = x(i)-h;
fh = feval(fcn, x, varargin{:});
fh = penalty_objective_function(x, fcn, penalty, varargin{:});
g0 = (f0-fh)/h;
end
if abs(g0)< 1e15

6
matlab/optimization/numgrad3.m
View File

@ -1,4 +1,4 @@
function [g, badg] = numgrad3(fcn,f0,x,epsilon,varargin)
function [g, badg] = numgrad3(fcn,f0,x,penalty,epsilon,varargin)
% Computes the gradient of the objective function fcn using a three points
% formula if possible.
%
@ -38,9 +38,9 @@ badg=0;
for i=1:n
xiold = x(i);
x(i) = xiold+h;
f1 = feval(fcn, x, varargin{:});
f1 = penalty_objective_function(x, fcn, penatly, varargin{:});
x(i) = xiold-h;
f2 = feval(fcn, x, varargin{:});
f2 = penalty_objective_function(x, fcn, penatly, varargin{:});
g0 = (f1-f2)/H;
if abs(g0)< 1e15
g(i)=g0;

6
matlab/optimization/numgrad3_.m
View File

@ -1,4 +1,4 @@
function [g, badg] = numgrad3_(fcn,f0,x,epsilon,varargin)
function [g, badg] = numgrad3_(fcn,f0,x,penalty,epsilon,varargin)
% Computes the gradient of the objective function fcn using a three points
% formula if possible.
%
@ -46,12 +46,12 @@ for i=1:n
xiold = x(i);
h = step_length_correction(xiold,scale,i)*delta;
x(i) = xiold + h;
[f1,junk1,junk2,cost_flag1] = feval(fcn, x, varargin{:});
[f1,junk1,junk2,cost_flag1] = penalty_objective_function(x, fcn, penalty, varargin{:});
if ~cost_flag1
fprintf('Gradient w.r.t. parameter number %3d (x=%16.8f,+h=%16.8f,f0=%16.8f,f1=%16.8f,f2=%16.8f,g0=%16.8f): penalty on the right!\n',i,xiold,h,f0,f1,f2,(f1 - f2) / (2*h))
end
x(i) = xiold - h;
[f2,junk2,junk2,cost_flag2] = feval(fcn, x, varargin{:});
[f2,junk2,junk2,cost_flag2] = penalty_objective_function(x, fcn, penalty, varargin{:});
if ~cost_flag2
fprintf('Gradient w.r.t. parameter number %3d (x=%16.8f,+h=%16.8f,f0=%16.8f,f1=%16.8f,f2=%16.8f,g0=%16.8f): penalty on the left!\n',i,xiold,h,f0,f1,f2,(f1 - f2) / (2*h))
end

10
matlab/optimization/numgrad5.m
View File

@ -1,4 +1,4 @@
function [g, badg] = numgrad5(fcn,f0,x,epsilon,varargin)
function [g, badg] = numgrad5(fcn,f0,x,penalty,epsilon,varargin)
% Computes the gradient of the objective function fcn using a five points
% formula if possible.
%
@ -40,13 +40,13 @@ badg = 0;
for i=1:n
xiold = x(i);
x(i) = xiold+h;
f1 = feval(fcn, x, varargin{:});
f1 = penalty_objective_function(x, fcn, penalty, varargin{:});
x(i) = xiold-h;
f2 = feval(fcn, x, varargin{:});
f2 = penalty_objective_function(x, fcn, penalty, varargin{:});
x(i) = xiold+2*h;
f3 = feval(fcn, x, varargin{:});
f3 = penalty_objective_function(x, fcn, penalty, varargin{:});
x(i) = xiold-2*h;
f4 = feval(fcn, x, varargin{:});
f4 = penalty_objective_function(x, fcn, penalty, varargin{:});
g0 = (8*(f1-f2)+f4-f3)/H;
if abs(g0)< 1e15
g(i) = g0;

8
matlab/optimization/numgrad5_.m
View File

@ -48,13 +48,13 @@ for i=1:n
xiold = x(i);
h = step_length_correction(xiold,scale,i)*delta;
x(i) = xiold+h;
[f1,junk1,junk2,cost_flag1] = feval(fcn, x, varargin{:});
[f1,junk1,junk2,cost_flag1] = feval(fcn, x, penalty, varargin{:});
x(i) = xiold-h;
[f2,junk1,junk2,cost_flag2] = feval(fcn, x, varargin{:});
[f2,junk1,junk2,cost_flag2] = feval(fcn, x, penalty, varargin{:});
x(i) = xiold+2*h;
[f3,junk1,junk2,cost_flag3] = feval(fcn, x, varargin{:});
[f3,junk1,junk2,cost_flag3] = feval(fcn, x, penalty, varargin{:});
x(i) = xiold-2*h;
[f4,junk1,junk2,cost_flag4] = feval(fcn, x, varargin{:});
[f4,junk1,junk2,cost_flag4] = feval(fcn, x, penalty, varargin{:});
if f0<f1 && f1<f3 && f0<f2 && f2<f4
g0 = 0;
else