From 837c812c1700e7abef7bafdfb22c79568b3073ce Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Hermes=29?= Date: Fri, 26 Feb 2016 22:06:49 +0100 Subject: [PATCH] Efficiency changes. + Use scalar product when possible, ie replace sum(x.*x) by x'*x, + Add parenthesis around scalar products, ie replace a*x'*x by a*(x'*x), + Removed global options_ in lnsrch1 (only used for solve_tolx). Note that the scalar products could be replaced by a loop for small problems. If the number of unknow is smaller than 70 the loop is faster (10 times faster if the number of unknowns is less than 10). --- matlab/lnsrch1.m | 24 ++++++++---------------- matlab/solve1.m | 6 +++--- matlab/solve_one_boundary.m | 6 +++--- matlab/solve_two_boundaries.m | 4 ++-- 4 files changed, 16 insertions(+), 24 deletions(-) diff --git a/matlab/lnsrch1.m b/matlab/lnsrch1.m index 80af34462..dbcbc5918 100644 --- a/matlab/lnsrch1.m +++ b/matlab/lnsrch1.m @@ -1,5 +1,5 @@ -function [x,f,fvec,check]=lnsrch1(xold,fold,g,p,stpmax,func,j1,j2,varargin) -% function [x,f,fvec,check]=lnsrch1(xold,fold,g,p,stpmax,func,j1,j2,varargin) +function [x,f,fvec,check]=lnsrch1(xold, fold, g, p, stpmax, func, j1, j2, tolx, varargin) +% function [x,f,fvec,check]=lnsrch1(xold,fold,g,p,stpmax,func,j1,j2,tolx,varargin) % Computes the optimal step by minimizing the residual sum of squares % % INPUTS @@ -11,6 +11,7 @@ function [x,f,fvec,check]=lnsrch1(xold,fold,g,p,stpmax,func,j1,j2,varargin) % func: name of the function % j1: equations index to be solved % j2: unknowns index +% tolx: tolerance parameter % varargin: list of arguments following j2 % % OUTPUTS @@ -23,7 +24,7 @@ function [x,f,fvec,check]=lnsrch1(xold,fold,g,p,stpmax,func,j1,j2,varargin) % SPECIAL REQUIREMENTS % none -% Copyright (C) 2001-2010 Dynare Team +% Copyright (C) 2001-2016 Dynare Team % % This file is part of Dynare. % @@ -40,15 +41,13 @@ function [x,f,fvec,check]=lnsrch1(xold,fold,g,p,stpmax,func,j1,j2,varargin) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -global options_ - alf = 1e-4 ; -tolx = options_.solve_tolx; alam = 1; x = xold; nn = length(j2); -summ = sqrt(sum(p.*p)) ; +summ = sqrt(p'*p); + if ~isfinite(summ) eq_number_string=[]; for ii=1:length(j1)-1 @@ -73,7 +72,7 @@ if ~isfinite(summ) end if summ > stpmax - p=p.*stpmax/summ ; + p = p*stpmax/summ ; end slope = g'*p ; @@ -90,19 +89,16 @@ while 1 check = 1 ; return end - x(j2) = xold(j2) + (alam*p) ; fvec = feval(func,x,varargin{:}) ; fvec = fvec(j1); - f = 0.5*fvec'*fvec ; - + f = 0.5*(fvec'*fvec) ; if any(isnan(fvec)) alam = alam/2 ; alam2 = alam ; f2 = f ; fold2 = fold ; else - if f <= fold+alf*alam*slope check = 0; break ; @@ -124,15 +120,11 @@ while 1 else tmplam = (-b+sqrt(disc))/(3*a) ; end - end - if tmplam > 0.5*alam tmplam = 0.5*alam; end - end - alam2 = alam ; f2 = f ; fold2 = fold ; diff --git a/matlab/solve1.m b/matlab/solve1.m index 5fecdd356..937909d09 100644 --- a/matlab/solve1.m +++ b/matlab/solve1.m @@ -23,7 +23,7 @@ function [x,check] = solve1(func,x,j1,j2,jacobian_flag,gstep,tolf,tolx,maxit,deb % SPECIAL REQUIREMENTS % none -% Copyright (C) 2001-2014 Dynare Team +% Copyright (C) 2001-2016 Dynare Team % % This file is part of Dynare. % @@ -62,7 +62,7 @@ if ~isempty(i) return; end -f = 0.5*fvec'*fvec ; +f = 0.5*(fvec'*fvec) ; if max(abs(fvec)) < tolf return ; @@ -108,7 +108,7 @@ for its = 1:maxit xold = x ; fold = f ; - [x,f,fvec,check]=lnsrch1(xold,fold,g,p,stpmax,func,j1,j2,varargin{:}); + [x, f, fvec, check] = lnsrch1(xold, fold, g, p, stpmax, func, j1, j2, tolx, varargin{:}); if debug disp([its f]) diff --git a/matlab/solve_one_boundary.m b/matlab/solve_one_boundary.m index b1a7ed793..2ce90c7bc 100644 --- a/matlab/solve_one_boundary.m +++ b/matlab/solve_one_boundary.m @@ -55,7 +55,7 @@ function [y, info] = solve_one_boundary(fname, y, x, params, steady_state, ... % none. % -% Copyright (C) 1996-2015 Dynare Team +% Copyright (C) 1996-2016 Dynare Team % % This file is part of Dynare. % @@ -247,10 +247,10 @@ for it_=start:incr:finish if is_dynamic [ya,f,r,check]=lnsrch1(y(it_,:)',f,g,p,stpmax, ... 'lnsrch1_wrapper_one_boundary',nn, ... - y_index_eq, y_index_eq, fname, y, x, params, steady_state, it_); + y_index_eq, options.solve_tolx, y_index_eq, fname, y, x, params, steady_state, it_); dx = ya' - y(it_, :); else - [ya,f,r,check]=lnsrch1(y,f,g,p,stpmax,fname,nn,y_index_eq,x, ... + [ya,f,r,check]=lnsrch1(y,f,g,p,stpmax,fname,nn,y_index_eq, options.solve_tolx, x, ... params, steady_state,0); dx = ya - y(y_index_eq); end diff --git a/matlab/solve_two_boundaries.m b/matlab/solve_two_boundaries.m index cf00fdacb..136d3e1ac 100644 --- a/matlab/solve_two_boundaries.m +++ b/matlab/solve_two_boundaries.m @@ -46,7 +46,7 @@ function [y, oo]= solve_two_boundaries(fname, y, x, params, steady_state, y_inde % none. % -% Copyright (C) 1996-2015 Dynare Team +% Copyright (C) 1996-2016 Dynare Team % % This file is part of Dynare. % @@ -307,7 +307,7 @@ while ~(cvg==1 || iter>maxit_) g = (ra'*g1a)'; f = 0.5*ra'*ra; p = -g1a\ra; - [yn,f,ra,check]=lnsrch1(ya,f,g,p,stpmax,'lnsrch1_wrapper_two_boundaries',nn,nn, fname, y, y_index,x, params, steady_state, periods, y_kmin, Blck_size,options.periods); + [yn,f,ra,check]=lnsrch1(ya,f,g,p,stpmax,'lnsrch1_wrapper_two_boundaries',nn,nn, options.solve_tolx, fname, y, y_index,x, params, steady_state, periods, y_kmin, Blck_size,options.periods); dx = ya - yn; y(1+y_kmin:periods+y_kmin,y_index)=reshape(yn',length(y_index),periods)'; end