Commit 6d0c4133 authored by Marco Ratto's avatar Marco Ratto
Browse files

added new robustness features to newrat. manually fixed conflicts

(cherry picked from commit 09555f4b3461cfd31bfc411ee222aadb0c606707)
parent e72e6fb3
......@@ -288,6 +288,8 @@ switch minimizer_algorithm
newratflag = options_.newrat.hess; %default
end
nit=options_.newrat.maxiter;
epsilon = options_.gradient_epsilon;
crit_rel = NaN;
Verbose = options_.newrat.verbosity;
Save_files = options_.newrat.Save_files;
if ~isempty(options_.optim_opt)
......@@ -303,8 +305,12 @@ switch minimizer_algorithm
else
newratflag=flag;
end
case 'NumgradEpsilon'
epsilon = options_list{i,2};
case 'TolFun'
crit = options_list{i,2};
case 'TolFunRel'
crit_rel = options_list{i,2};
case 'verbosity'
Verbose = options_list{i,2};
case 'SaveFiles'
......@@ -318,9 +324,10 @@ switch minimizer_algorithm
Save_files = 0;
Verbose = 0;
end
hess_info.crit_rel=crit_rel;
hess_info.gstep=options_.gstep;
hess_info.htol = 1.e-4;
hess_info.h1=options_.gradient_epsilon*ones(n_params,1);
hess_info.h1=epsilon*ones(n_params,1);
% here we force 7th input argument (flagg) to be 0, since outer product
% gradient Hessian is handled in dynare_estimation_1
[opt_par_values,hessian_mat,gg,fval,invhess,new_rat_hess_info] = newrat(objective_function,start_par_value,bounds,analytic_grad,crit,nit,0,Verbose,Save_files,hess_info,prior_information.p2,options_.gradient_epsilon,parameter_names,varargin{:}); %hessian_mat is the plain outer product gradient Hessian
......
......@@ -72,15 +72,77 @@ 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(i)
[f0, x, ~, ~] = csminit1(func0,x,penalty,f0,gg,0,diag(hh),Verbose,varargin{:});
[ff, xx,fcount,retcode] = csminit1(func0,x,penalty,f0,gg,0,diag(hh),Verbose,varargin{:});
if retcode
% do_nothing=true;
xa=transpose(linspace(x(i)/2, x(i)*3/2, 7));
for k=1:7
xh1(i)=xa(k);
fa(k,1) = penalty_objective_function(xh1,func0,penalty,varargin{:});
end
b=[ones(7,1) xa xa.*xa./2]\fa;
gg(i)=x(i)*b(3)+b(2);
hh(i)=1/b(3);
[ff2, xx2, fcount2, retcode2] = csminit1(func0,x,penalty,f0,gg,0,diag(hh),Verbose,varargin{:});
if ff2<ff
ff=ff2;
xx=xx2;
end
if min(fa)<ff
[ff, im]=min(fa);
xx(i)=xa(im);
end
end
ig(i)=1;
if not(isequal(xx , check_bounds(xx,bounds)))
xx = check_bounds(xx,bounds);
if xx(i)<x(i)
% lower bound
xx(i) = min(xx(i)+h1(i), 0.5*(xx(i)+x(i)));
else
% upper bound
xx(i) = max(xx(i)-h1(i), 0.5*(xx(i)+x(i)));
end
[ff,exit_flag]=penalty_objective_function(xx,func0,penalty,varargin{:});
if exit_flag~=1
disp('last step exited with bad status!')
elseif ff<f0
f0=ff;
x=xx;
end
else
% check improvement wrt predicted one
if abs(f0-ff) < abs(gg(i)*(hh(i)*gg(i))/2/100) || abs(x(i)-xx(i))<1.e-10
[ff1, xx1, fcount, retcode] = csminit1(func0,x,penalty,f0,-gg,0,diag(hh),Verbose,varargin{:});
if not(isequal(xx1 , check_bounds(xx1,bounds)))
xx1 = check_bounds(xx1,bounds);
if xx1(i)<x(i)
% lower bound
xx1(i) = min(xx1(i)+h1(i), 0.5*(xx1(i)+x(i)));
else
% upper bound
xx1(i) = max(xx1(i)-h1(i), 0.5*(xx1(i)+x(i)));
end
[ff1,exit_flag]=penalty_objective_function(xx1,func0,penalty,varargin{:});
if exit_flag~=1
disp('last step exited with bad status!')
end
end
if ff1<ff
ff=ff1;
xx=xx1;
end
end
f0=ff;
x=xx;
end
if Verbose
fprintf(['Done for param %s = %8.4f\n'],parameter_names{i},x(i))
fprintf(['Done for param %s = %8.4f; f = %8.4f\n'],parameter_names{i},x(i),f0)
end
end
xh1=x;
end
x = check_bounds(x,bounds);
% penalty=f0;
if Save_files
save('gstep.mat','x','h1','f0')
end
......@@ -96,10 +158,10 @@ function x = check_bounds(x,bounds)
inx = find(x>=bounds(:,2));
if ~isempty(inx)
x(inx) = bounds(inx,2)-eps;
x(inx) = bounds(inx,2)-1.e-10;
end
inx = find(x<=bounds(:,1));
if ~isempty(inx)
x(inx) = bounds(inx,1)+eps;
x(inx) = bounds(inx,1)+1.e-10;
end
......@@ -78,11 +78,16 @@ else
end
hess_info.h1 = min(hess_info.h1,0.5.*hmax);
hess_info.h1 = min(hess_info.h1,0.9.*hmax);
if htol0<hess_info.htol
hess_info.htol=htol0;
end
if not(isnan(hess_info.crit_rel))
htol1=hess_info.htol;
hess_info.htol=abs(hess_info.crit_rel*f0);
end
xh1=x;
f1=zeros(size(f0,1),n);
f_1=f1;
......@@ -106,22 +111,39 @@ while i<n
ic=0;
icount = 0;
h0=hess_info.h1(i);
while (abs(dx(it))<0.5*hess_info.htol || abs(dx(it))>(3*hess_info.htol)) && icount<10 && ic==0
istoobig=false;
istoosmall=false;
give_it_up=false;
while (abs(dx(it))<0.5*hess_info.htol || abs(dx(it))>(3*hess_info.htol)) && icount<10 && ic==0 && hmax(i)>2*1.e-10 && not(give_it_up)
icount=icount+1;
istoobig(it)=false;
istoosmall(it)=false;
if abs(dx(it))<0.5*hess_info.htol
istoosmall(it)=true;
if hess_info.h1(i)==0.9*hmax(i)% || hess_info.h1(i)==0.3*abs(x(i))
give_it_up=true;
end
if abs(dx(it)) ~= 0
hess_info.h1(i)=min(max(1.e-10,0.3*abs(x(i))), 0.9*hess_info.htol/abs(dx(it))*hess_info.h1(i));
% htmp=min(max(1.e-10,0.3*abs(x(i))), 0.9*hess_info.htol/abs(dx(it))*hess_info.h1(i));
htmp=0.9*hess_info.htol/abs(dx(it))*hess_info.h1(i);
else
hess_info.h1(i)=2.1*hess_info.h1(i);
htmp=2.1*hess_info.h1(i);
end
hess_info.h1(i) = min(hess_info.h1(i),0.5*hmax(i));
hess_info.h1(i) = max(hess_info.h1(i),1.e-10);
htmp = min(htmp,0.9*hmax(i));
if any(h0(istoobig(1:it))) %&& htmp>=min(h0(istoobig(1:it)))
htmp = 0.5*(min(h0(istoobig))+max(h0(istoosmall)));
end
hess_info.h1(i) = max(htmp,1.e-10);
xh1(i)=x(i)+hess_info.h1(i);
[fx,exit_flag,ffx]=penalty_objective_function(xh1,func,penalty,varargin{:});
end
if abs(dx(it))>(3*hess_info.htol)
hess_info.h1(i)= hess_info.htol/abs(dx(it))*hess_info.h1(i);
hess_info.h1(i) = max(hess_info.h1(i),1e-10);
istoobig(it)=true;
htmp= hess_info.htol/abs(dx(it))*hess_info.h1(i);
if any(h0(istoosmall(1:it))) %&& htmp<=max(h0(istoosmall(1:it)))
htmp = 0.5*(min(h0(istoobig))+max(h0(istoosmall)));
end
hess_info.h1(i) = max(htmp,1e-10);
xh1(i)=x(i)+hess_info.h1(i);
[fx,exit_flag,ffx]=penalty_objective_function(xh1,func,penalty,varargin{:});
iter=0;
......@@ -136,11 +158,27 @@ while i<n
it=it+1;
dx(it)=(fx-f0);
h0(it)=hess_info.h1(i);
if (hess_info.h1(i)<1.e-12*min(1,h2(i)) && hess_info.h1(i)<0.5*hmax(i))
if (hess_info.h1(i)<1.e-12*min(1,h2(i)) && hess_info.h1(i)<0.9*hmax(i))
ic=1;
hcheck=1;
end
end
if icount == 10 || hess_info.h1(i)==1.e-10
istoobig(it)=false;
istoosmall(it)=false;
if abs(dx(it))<0.5*hess_info.htol
istoosmall(it)=true;
end
if abs(dx(it))>(3*hess_info.htol)
istoobig(it)=true;
end
if any(istoobig) && (istoosmall(it) || istoobig(it))
% always better to be wrong from above, to avoid numerical noise
[ddx, ij]=min(dx(istoobig));
fx=f0+ddx;
hess_info.h1(i)=h0(ij);
end
end
f1(:,i)=fx;
if outer_product_gradient
if any(isnan(ffx)) || isempty(ffx)
......@@ -265,4 +303,6 @@ else
hh1 = [];
end
htol1=hhtol;
if isnan(hess_info.crit_rel)
htol1=hhtol;
end
\ No newline at end of file
......@@ -94,6 +94,7 @@ fval=fval0;
outer_product_gradient=1;
if isempty(hh)
penalty=fval0;
[dum, gg, htol0, igg, hhg, h1, hess_info]=mr_hessian(x,func0,penalty,flagit,htol,hess_info,bounds,prior_std,Save_files,varargin{:});
if isempty(dum)
outer_product_gradient=0;
......@@ -176,9 +177,25 @@ while norm(gg)>gtol && check==0 && jit<nit
iggx(ig_pos,ig_pos) = inv( hhx(ig_pos,ig_pos) );
[fvala,x0,fc,retcode] = csminit1(func0,x0,penalty,fval,ggx,0,iggx,Verbose,varargin{:});
end
x0 = check_bounds(x0,bounds);
if not(isequal(x0 , check_bounds(x0,bounds)))
x0 = check_bounds(x0,bounds);
[fvala,exit_flag]=penalty_objective_function(x0,func0,penalty,varargin{:});
if exit_flag==1
penalty=fvala;
else
disp('last step exited with bad status!')
end
end
[fvala, x0, ig] = mr_gstep(h1,x0,bounds,func0,penalty,htol0,Verbose,Save_files,gradient_epsilon, parameter_names,varargin{:});
x0 = check_bounds(x0,bounds);
if not(isequal(x0 , check_bounds(x0,bounds)))
x0 = check_bounds(x0,bounds);
[fvala,exit_flag]=penalty_objective_function(x0,func0,penalty,varargin{:});
if exit_flag==1
penalty=fvala;
else
disp('last step exited with bad status!')
end
end
if Save_files
nig=[nig ig];
end
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment