clean-up of commented lines in optimizer number 5;
eliminated globals and persistent from mr_gstep;time-shift
Marco Ratto
parent
da664f5a26
commit
aac282d371
|
|
@ -1,5 +1,5 @@
|
|||
function [f0, x, ig] = mr_gstep(init,x,func0,htol0,varargin)
|
||||
% function [f0, x] = mr_gstep(init,x,func0,htol0,varargin)
|
||||
function [f0, x, ig] = mr_gstep(h1,x,func0,htol0,varargin)
|
||||
% function [f0, x, ig] = mr_gstep(h1,x,func0,htol0,varargin)
|
||||
%
|
||||
% Gibbs type step in optimisation
|
||||
|
||||
|
|
@ -20,17 +20,8 @@ function [f0, x, ig] = mr_gstep(init,x,func0,htol0,varargin)
|
|||
% You should have received a copy of the GNU General Public License
|
||||
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
global bayestopt_ options_
|
||||
persistent h1
|
||||
|
||||
n=size(x,1);
|
||||
|
||||
if init,
|
||||
gstep_ = options_.gstep;
|
||||
% h1=max(abs(x),sqrt(gstep_)*ones(n,1))*eps^(1/4);
|
||||
h1=options_.gradient_epsilon*ones(n,1);
|
||||
return
|
||||
end
|
||||
if nargin<4,
|
||||
htol = 1.e-6;
|
||||
else
|
||||
|
|
@ -38,12 +29,11 @@ else
|
|||
end
|
||||
func = str2func(func0);
|
||||
f0=feval(func,x,varargin{:});
|
||||
h2=bayestopt_.ub-bayestopt_.lb;
|
||||
|
||||
xh1=x;
|
||||
f1=zeros(size(f0,1),n);
|
||||
f_1=f1;
|
||||
%for i=1:n,
|
||||
|
||||
i=0;
|
||||
ig=zeros(n,1);
|
||||
while i<n,
|
||||
|
|
@ -53,78 +43,12 @@ while i<n,
|
|||
dx=[];
|
||||
xh1(i)=x(i)+h1(i);
|
||||
fx = feval(func,xh1,varargin{:});
|
||||
it=1;
|
||||
dx=(fx-f0);
|
||||
ic=0;
|
||||
% if abs(dx)>(2*htol),
|
||||
% c=mr_nlincon(xh1,varargin{:});
|
||||
% while c
|
||||
% h1(i)=h1(i)*0.9;
|
||||
% xh1(i)=x(i)+h1(i);
|
||||
% c=mr_nlincon(xh1,varargin{:});
|
||||
% ic=1;
|
||||
% end
|
||||
% if ic,
|
||||
% fx = feval(func,xh1,varargin{:});
|
||||
% dx=(fx-f0);
|
||||
% end
|
||||
% end
|
||||
|
||||
icount = 0;
|
||||
h0=h1(i);
|
||||
% while (abs(dx(it))<0.5*htol || abs(dx(it))>(2*htol)) && icount<10 && ic==0,
|
||||
while (abs(dx(it))<0.5*htol) && icount<10 && ic==0,
|
||||
%while abs(dx(it))<0.5*htol && icount< 10 && ic==0,
|
||||
icount=icount+1;
|
||||
if abs(dx(it)) ~= 0,
|
||||
if abs(dx(it))<0.5*htol
|
||||
h1(i)=min(0.3*abs(x(i)), 0.9*htol/abs(dx(it))*h1(i));
|
||||
xh1(i)=x(i)+h1(i);
|
||||
% c=mr_nlincon(xh1,varargin{:});
|
||||
% while c
|
||||
% h1(i)=h1(i)*0.9;
|
||||
% xh1(i)=x(i)+h1(i);
|
||||
% c=mr_nlincon(xh1,varargin{:});
|
||||
% ic=1;
|
||||
% end
|
||||
end
|
||||
% if abs(dx(it))>(2*htol),
|
||||
% h1(i)= htol/abs(dx(it))*h1(i);
|
||||
% xh1(i)=x(i)+h1(i);
|
||||
% end
|
||||
try
|
||||
fx = feval(func,xh1,varargin{:});
|
||||
catch
|
||||
fx=1.e8;
|
||||
end
|
||||
it=it+1;
|
||||
dx(it)=(fx-f0);
|
||||
h0(it)=h1(i);
|
||||
if h1(i)<1.e-12*min(1,h2(i)),
|
||||
ic=1;
|
||||
hcheck=1;
|
||||
end
|
||||
else
|
||||
h1(i)=1;
|
||||
ic=1;
|
||||
end
|
||||
end
|
||||
f1(:,i)=fx;
|
||||
xh1(i)=x(i)-h1(i);
|
||||
% c=mr_nlincon(xh1,varargin{:});
|
||||
% ic=0;
|
||||
% while c
|
||||
% h1(i)=h1(i)*0.9;
|
||||
% xh1(i)=x(i)-h1(i);
|
||||
% c=mr_nlincon(xh1,varargin{:});
|
||||
% ic = 1;
|
||||
% end
|
||||
|
||||
fx = feval(func,xh1,varargin{:});
|
||||
f_1(:,i)=fx;
|
||||
% if ic,
|
||||
% xh1(i)=x(i)+h1(i);
|
||||
% f1(:,i)=feval(func,xh1,varargin{:});
|
||||
% end
|
||||
|
||||
if hcheck && htol<1,
|
||||
htol=min(1,max(min(abs(dx))*2,htol*10));
|
||||
h1(i)=h10;
|
||||
|
|
@ -134,11 +58,12 @@ while i<n,
|
|||
gg=zeros(size(x));
|
||||
hh=gg;
|
||||
gg(i)=(f1(i)'-f_1(i)')./(2.*h1(i));
|
||||
if abs(f1(i)+f_1(i)-2*f0)>1.e-12,
|
||||
hh(i) = abs(1/( (f1(i)+f_1(i)-2*f0)./(h1(i)*h1(i)) ));
|
||||
else
|
||||
hh(i) = 1;
|
||||
end
|
||||
hh(i) = 1/max(1.e-9,abs( (f1(i)+f_1(i)-2*f0)./(h1(i)*h1(i)) ));
|
||||
% if abs(f1(i)+f_1(i)-2*f0)>1.e-12,
|
||||
% hh(i) = abs(1/( (f1(i)+f_1(i)-2*f0)./(h1(i)*h1(i)) ));
|
||||
% else
|
||||
% hh(i) = 1;
|
||||
% end
|
||||
|
||||
if gg(i)*(hh(i)*gg(i))/2 > htol,
|
||||
[f0 x fc retcode] = csminit(func0,x,f0,gg,0,diag(hh),varargin{:});
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
function [hessian_mat, gg, htol1, ihh, hh_mat0] = mr_hessian(init,x,func,hflag,htol0,varargin)
|
||||
% [hessian_mat, gg, htol1, ihh, hh_mat0] = mr_hessian(init,x,func,hflag,htol0,varargin)
|
||||
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)
|
||||
%
|
||||
% numerical gradient and Hessian, with 'automatic' check of numerical
|
||||
% error
|
||||
|
|
@ -47,10 +47,8 @@ persistent h1 htol
|
|||
n=size(x,1);
|
||||
if init,
|
||||
gstep_=options_.gstep;
|
||||
htol = 1.e-4;
|
||||
%h1=max(abs(x),gstep_*ones(n,1))*eps^(1/3);
|
||||
%h1=max(abs(x),sqrt(gstep_)*ones(n,1))*eps^(1/6);
|
||||
% h1=max(abs(x),sqrt(gstep_)*ones(n,1))*eps^(1/4);
|
||||
htol = 1.e-4;
|
||||
%h1=max(abs(x),sqrt(gstep_)*ones(n,1))*eps^(1/4);
|
||||
h1=options_.gradient_epsilon*ones(n,1);
|
||||
return,
|
||||
end
|
||||
|
|
@ -62,7 +60,7 @@ hmax=min(hmax,x-bayestopt_.lb);
|
|||
|
||||
h1 = min(h1,0.5.*hmax);
|
||||
|
||||
if htol0<htol,
|
||||
if htol0<htol,
|
||||
htol=htol0;
|
||||
end
|
||||
xh1=x;
|
||||
|
|
@ -70,14 +68,13 @@ f1=zeros(size(f0,1),n);
|
|||
f_1=f1;
|
||||
ff1=zeros(size(ff0));
|
||||
ff_1=ff1;
|
||||
ggh=zeros(size(ff0,1),n);
|
||||
|
||||
%for i=1:n,
|
||||
i=0;
|
||||
while i<n,
|
||||
i=i+1;
|
||||
h10=h1(i);
|
||||
hcheck=0;
|
||||
dx=[];
|
||||
xh1(i)=x(i)+h1(i);
|
||||
try
|
||||
[fx, ffx]=feval(func,xh1,varargin{:});
|
||||
|
|
@ -87,27 +84,12 @@ while i<n,
|
|||
it=1;
|
||||
dx=(fx-f0);
|
||||
ic=0;
|
||||
% if abs(dx)>(2*htol),
|
||||
% c=mr_nlincon(xh1,varargin{:});
|
||||
% while c
|
||||
% h1(i)=h1(i)*0.9;
|
||||
% xh1(i)=x(i)+h1(i);
|
||||
% c=mr_nlincon(xh1,varargin{:});
|
||||
% ic=1;
|
||||
% end
|
||||
% if ic,
|
||||
% [fx, ffx]=feval(func,xh1,varargin{:});
|
||||
% dx=(fx-f0);
|
||||
% end
|
||||
% end
|
||||
|
||||
icount = 0;
|
||||
h0=h1(i);
|
||||
% while (abs(dx(it))<0.5*htol || abs(dx(it))>(2*htol)) && icount<10 && ic==0,
|
||||
while (abs(dx(it))<0.5*htol) && icount<10 && ic==0,
|
||||
while (abs(dx(it))<0.5*htol || abs(dx(it))>(3*htol)) && icount<10 && ic==0,
|
||||
%while abs(dx(it))<0.5*htol && icount< 10 && ic==0,
|
||||
icount=icount+1;
|
||||
%if abs(dx(it)) ~= 0,
|
||||
if abs(dx(it))<0.5*htol
|
||||
if abs(dx(it)) ~= 0,
|
||||
h1(i)=min(max(1.e-10,0.3*abs(x(i))), 0.9*htol/abs(dx(it))*h1(i));
|
||||
|
|
@ -116,88 +98,50 @@ while i<n,
|
|||
end
|
||||
h1(i) = min(h1(i),0.5*hmax(i));
|
||||
xh1(i)=x(i)+h1(i);
|
||||
% c=mr_nlincon(xh1,varargin{:});
|
||||
% while c
|
||||
% h1(i)=h1(i)*0.9;
|
||||
% xh1(i)=x(i)+h1(i);
|
||||
% c=mr_nlincon(xh1,varargin{:});
|
||||
% ic=1;
|
||||
% end
|
||||
try
|
||||
try
|
||||
[fx, ffx]=feval(func,xh1,varargin{:});
|
||||
catch
|
||||
fx=1.e8;
|
||||
end
|
||||
end
|
||||
% if abs(dx(it))>(2*htol),
|
||||
% h1(i)= htol/abs(dx(it))*h1(i);
|
||||
% xh1(i)=x(i)+h1(i);
|
||||
% try
|
||||
% [fx, ffx]=feval(func,xh1,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{:});
|
||||
% ic=1;
|
||||
% end
|
||||
% end
|
||||
if abs(dx(it))>(3*htol),
|
||||
h1(i)= htol/abs(dx(it))*h1(i);
|
||||
xh1(i)=x(i)+h1(i);
|
||||
try
|
||||
[fx, ffx]=feval(func,xh1,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{:});
|
||||
ic=1;
|
||||
end
|
||||
end
|
||||
it=it+1;
|
||||
dx(it)=(fx-f0);
|
||||
h0(it)=h1(i);
|
||||
if h1(i)<1.e-12*min(1,h2(i)) && h1(i)<0.5*hmax(i),
|
||||
if (h1(i)<1.e-12*min(1,h2(i)) && h1(i)<0.5*hmax(i)),% || (icount==10 && abs(dx(it))>(3*htol)),
|
||||
ic=1;
|
||||
hcheck=1;
|
||||
end
|
||||
%else
|
||||
% h1(i)=1;
|
||||
% ic=1;
|
||||
%end
|
||||
end
|
||||
% if (it>2 && dx(1)<10^(log10(htol)/2)) ,
|
||||
% [dum, is]=sort(h0);
|
||||
% if find(diff(sign(diff(dx(is)))));
|
||||
% hcheck=1;
|
||||
% end
|
||||
% elseif (it>3 && dx(1)>10^(log10(htol)/2)) ,
|
||||
% [dum, is]=sort(h0);
|
||||
% if find(diff(sign(diff(dx(is(1:end-1))))));
|
||||
% hcheck=1;
|
||||
% end
|
||||
% end
|
||||
f1(:,i)=fx;
|
||||
if any(isnan(ffx)),
|
||||
ff1=ones(size(ff0)).*fx/length(ff0);
|
||||
else
|
||||
ff1=ffx;
|
||||
end
|
||||
% if hflag, % two point based derivatives
|
||||
xh1(i)=x(i)-h1(i);
|
||||
% c=mr_nlincon(xh1,varargin{:});
|
||||
% ic=0;
|
||||
% while c
|
||||
% h1(i)=h1(i)*0.9;
|
||||
% xh1(i)=x(i)-h1(i);
|
||||
% c=mr_nlincon(xh1,varargin{:});
|
||||
% ic = 1;
|
||||
% end
|
||||
[fx, ffx]=feval(func,xh1,varargin{:});
|
||||
f_1(:,i)=fx;
|
||||
if any(isnan(ffx)),
|
||||
ff_1=ones(size(ff0)).*fx/length(ff0);
|
||||
else
|
||||
ff_1=ffx;
|
||||
end
|
||||
% if ic,
|
||||
% xh1(i)=x(i)+h1(i);
|
||||
% [f1(:,i), ff1]=feval(func,xh1,varargin{:});
|
||||
% end
|
||||
ggh(:,i)=(ff1-ff_1)./(2.*h1(i));
|
||||
% else
|
||||
% ggh(:,i)=(ff1-ff0)./h1(i);
|
||||
% end
|
||||
xh1(i)=x(i)-h1(i);
|
||||
[fx, ffx]=feval(func,xh1,varargin{:});
|
||||
f_1(:,i)=fx;
|
||||
if any(isnan(ffx)),
|
||||
ff_1=ones(size(ff0)).*fx/length(ff0);
|
||||
else
|
||||
ff_1=ffx;
|
||||
end
|
||||
ggh(:,i)=(ff1-ff_1)./(2.*h1(i));
|
||||
xh1(i)=x(i);
|
||||
if hcheck && htol<1,
|
||||
htol=min(1,max(min(abs(dx))*2,htol*10));
|
||||
|
|
@ -211,11 +155,7 @@ h_1=h1;
|
|||
xh1=x;
|
||||
xh_1=xh1;
|
||||
|
||||
% if hflag,
|
||||
gg=(f1'-f_1')./(2.*h1);
|
||||
% else
|
||||
% gg=(f1'-f0)./h1;
|
||||
% end
|
||||
gg=(f1'-f_1')./(2.*h1);
|
||||
|
||||
if hflag==2,
|
||||
gg=(f1'-f_1')./(2.*h1);
|
||||
|
|
@ -232,29 +172,8 @@ if hflag==2,
|
|||
xh1(j)=x(j)+h_1(j);
|
||||
xh_1(i)=x(i)-h1(i);
|
||||
xh_1(j)=x(j)-h_1(j);
|
||||
%hessian_mat(:,(i-1)*n+j)=-(-feval(func,xh1,varargin{:})-feval(func,xh_1,varargin{:})+temp(:,i)+temp(:,j))./(2*h1(i)*h_1(j));
|
||||
%temp1 = feval(func,xh1,varargin{:});
|
||||
% c=mr_nlincon(xh1,varargin{:});
|
||||
% lam=1;
|
||||
% while c,
|
||||
% lam=lam*0.9;
|
||||
% xh1(i)=x(i)+h1(i)*lam;
|
||||
% xh1(j)=x(j)+h_1(j)*lam;
|
||||
% %disp( ['hessian warning cross ', num2str(c) ]),
|
||||
% c=mr_nlincon(xh1,varargin{:});
|
||||
% end
|
||||
% temp1 = f0+(feval(func,xh1,varargin{:})-f0)/lam;
|
||||
temp1 = feval(func,xh1,varargin{:});
|
||||
|
||||
% c=mr_nlincon(xh_1,varargin{:});
|
||||
% while c,
|
||||
% lam=lam*0.9;
|
||||
% xh_1(i)=x(i)-h1(i)*lam;
|
||||
% xh_1(j)=x(j)-h_1(j)*lam;
|
||||
% %disp( ['hessian warning cross ', num2str(c) ]),
|
||||
% c=mr_nlincon(xh_1,varargin{:});
|
||||
% end
|
||||
% temp2 = f0+(feval(func,xh_1,varargin{:})-f0)/lam;
|
||||
temp2 = feval(func,xh_1,varargin{:});
|
||||
|
||||
hessian_mat(:,(i-1)*n+j)=-(-temp1 -temp2+temp(:,i)+temp(:,j))./(2*h1(i)*h_1(j));
|
||||
|
|
@ -285,8 +204,8 @@ gga=ggh.*kron(ones(size(ff1)),2.*h1'); % re-scaled gradient
|
|||
hh_mat=gga'*gga; % rescaled outer product hessian
|
||||
hh_mat0=ggh'*ggh; % outer product hessian
|
||||
A=diag(2.*h1); % rescaling matrix
|
||||
igg=inv(hh_mat); % inverted rescaled outer product hessian
|
||||
ihh=A'*igg*A; % inverted outer product hessian
|
||||
% igg=inv(hh_mat); % inverted rescaled outer product hessian
|
||||
ihh=A'*(hh_mat\A); % inverted outer product hessian
|
||||
if hflag>0 && min(eig(reshape(hessian_mat,n,n)))>0,
|
||||
hh0 = A*reshape(hessian_mat,n,n)*A'; %rescaled second order derivatives
|
||||
hh = reshape(hessian_mat,n,n); %rescaled second order derivatives
|
||||
|
|
@ -294,7 +213,7 @@ if hflag>0 && min(eig(reshape(hessian_mat,n,n)))>0,
|
|||
sd=sqrt(diag(hh_mat)); %rescaled 'standard errors' using outer product
|
||||
hh_mat=hh_mat./(sd*sd').*(sd0*sd0'); %rescaled inverse outer product with 'true' std's
|
||||
igg=inv(hh_mat); % rescaled outer product hessian with 'true' std's
|
||||
ihh=A'*igg*A; % inverted outer product hessian
|
||||
ihh=A'*(hh_mat\A); % inverted outer product hessian
|
||||
hh_mat0=inv(A)'*hh_mat*inv(A); % outer product hessian with 'true' std's
|
||||
sd=sqrt(diag(ihh)); %standard errors
|
||||
sdh=sqrt(1./diag(hh)); %diagonal standard errors
|
||||
|
|
|
|||
|
|
@ -60,11 +60,11 @@ func = str2func(func0);
|
|||
fval0=feval(func,x,varargin{:});
|
||||
fval=fval0;
|
||||
% initialize mr_gstep and mr_hessian
|
||||
mr_gstep(1,x);
|
||||
% mr_gstep(1,x);
|
||||
mr_hessian(1,x);
|
||||
|
||||
if isempty(hh)
|
||||
[dum, gg, htol0, igg, hhg]=mr_hessian(0,x,func_hh,flagit,htol,varargin{:});
|
||||
[dum, gg, htol0, igg, hhg, h1]=mr_hessian(0,x,func_hh,flagit,htol,varargin{:});
|
||||
hh0 = reshape(dum,nx,nx);
|
||||
hh=hhg;
|
||||
if min(eig(hh0))<0,
|
||||
|
|
@ -126,17 +126,12 @@ while norm(gg)>gtol && check==0 && jit<nit,
|
|||
iggx(find(ig),find(ig)) = inv( hhx(find(ig),find(ig)) );
|
||||
[fvala x0 fc retcode] = csminit(func0,x0,fval,ggx,0,iggx,varargin{:});
|
||||
end
|
||||
[fvala, x0, ig] = mr_gstep(0,x0,func0,htol,varargin{:});
|
||||
[fvala, x0, ig] = mr_gstep(h1,x0,func0,htol,varargin{:});
|
||||
% if length(find(ig))==0,
|
||||
% [fvala, x0, ig] = mr_gstep(0,x0,func0,htol/10,varargin{:});
|
||||
% [fvala, x0, ig] = mr_gstep(h1,x0,func0,htol/10,varargin{:});
|
||||
% end
|
||||
nig=[nig ig];
|
||||
if (fval-fvala)<gibbstol*(fval0(icount)-fval),
|
||||
igibbs=0;
|
||||
disp('Last sequence of univariate step, gain too small!!')
|
||||
else
|
||||
disp('Sequence of univariate steps!!')
|
||||
end
|
||||
disp('Sequence of univariate steps!!')
|
||||
fval=fvala;
|
||||
end
|
||||
if (fval0(icount)-fval)<ftol && flagit==0,
|
||||
|
|
@ -163,16 +158,13 @@ while norm(gg)>gtol && check==0 && jit<nit,
|
|||
xparam1=x0;
|
||||
x(:,icount+1)=xparam1;
|
||||
fval0(icount+1)=fval;
|
||||
%if (fval0(icount)-fval)<ftol*ftol && flagg==1;,
|
||||
mr_gstep(1,x);
|
||||
mr_hessian(1,x);
|
||||
if (fval0(icount)-fval)<ftol,
|
||||
disp('No further improvement is possible!')
|
||||
check=1;
|
||||
if flagit==2,
|
||||
hh=hh0;
|
||||
elseif flagg>0,
|
||||
[dum, gg, htol0, igg, hhg]=mr_hessian(0,xparam1,func_hh,flagg,ftol0,varargin{:});
|
||||
[dum, gg, htol0, igg, hhg,h1]=mr_hessian(0,xparam1,func_hh,flagg,ftol0,varargin{:});
|
||||
if flagg==2,
|
||||
hh = reshape(dum,nx,nx);
|
||||
ee=eig(hh);
|
||||
|
|
@ -202,9 +194,10 @@ while norm(gg)>gtol && check==0 && jit<nit,
|
|||
disp(['Ftol ',num2str(ftol)])
|
||||
disp(['Htol ',num2str(htol0)])
|
||||
|
||||
if df<htol0,
|
||||
htol=max(htol_base,df/10);
|
||||
end
|
||||
% if df<htol0,
|
||||
% htol=max(htol_base,df/10);
|
||||
% end
|
||||
htol=htol_base;
|
||||
|
||||
if norm(x(:,icount)-xparam1)>1.e-12,
|
||||
try
|
||||
|
|
@ -212,7 +205,7 @@ while norm(gg)>gtol && check==0 && jit<nit,
|
|||
catch
|
||||
save m1.mat x fval0 nig
|
||||
end
|
||||
[dum, gg, htol0, igg, hhg]=mr_hessian(0,xparam1,func_hh,flagit,htol,varargin{:});
|
||||
[dum, gg, htol0, igg, hhg, h1]=mr_hessian(0,xparam1,func_hh,flagit,htol,varargin{:});
|
||||
if htol0>htol, %ftol,
|
||||
%ftol=htol0;
|
||||
htol=htol0;
|
||||
|
|
|
|||
Loading…
Reference in New Issue