Speed up vectorized code automatic expansion logical indexing
2 ビュー (過去 30 日間)
古いコメントを表示
Alessandro Maria Marco
2022 年 10 月 19 日
コメント済み: Alessandro Maria Marco
2022 年 10 月 20 日
I have to compute the maximum of a 2-dimensional array repeatedly in a nested loop and I am trying to speed up the code, also using automatic expansion.
Basically I have an array F(l,a',a,z) and for each a and z (these are the loop variables) I find the maximum over l and a' and store it in some arrays. Consider that to compute the maximand I call a function. I already use automatic expansion in some places (see comments in the MWE) and I use logical indexing inside the function.
I attach a minimum working example so that anyone can run the problem. I tested it with Matlab online in this forum and it runs. Any help is greatly appreciated!
clear
clc
close all
na = 300;
nz = 66;
nl = 50;
% Generate fake data
l_grid = linspace(0.001,0.999,nl)' ;
a_grid = linspace(1e-6,250,na)';
z_grid = linspace(0.02,25,nz)';
% pi_z is a transition matrix
pi_z = rand(nz,nz);
pi_z = pi_z./sum(pi_z,2);
V_next = rand(na,nz);
% Initialize output of the loop
V = zeros(na,nz);
pol_ap = zeros(na,nz);
pol_l = zeros(na,nz);
% Choice variables: (l,a')
l_choice = l_grid; %(nl,1)
aprime_choice = a_grid'; %(1,na)
%% Numerical computation starts here
tic
% Expected value EV(a') given z,theta and z' integrated out
EV = V_next*pi_z'; % V(a',z')*pi(z,z')T ==> EV(a',z)
% Loop over shock z today
for z_c = 1:nz
EVz = EV(:,z_c)';
% Loop over assets today
for a_c = 1:na
a_val = a_grid(a_c); % assets today, scalar
z_val = z_grid(z_c); % shock today, scalar
RHS_mat = myfun(l_choice,aprime_choice,a_val,z_val)+EVz;
[V(a_c,z_c),opt_ind] = max(RHS_mat(:));
[l_opt_ind,ap_opt_ind] = ind2sub([nl,na],opt_ind);
aprime_opt = a_grid(ap_opt_ind);
l_opt = l_grid(l_opt_ind);
pol_ap(a_c,z_c) = aprime_opt;
pol_l(a_c,z_c) = l_opt;
end % end loop over assets a
end % end loop over z
toc
function F = myfun(l,aprime,a_val,z_val)
aux = 3.5;
cons = a_val+l*z_val - aprime; % Use automatic expansion, dim: (nl,na)
indx = cons<=0; % NOT feasible!
F = log(cons)-l.^aux/aux; % Use automatic expansion
F(indx) = -inf; % Logical indexing
end
6 件のコメント
dpb
2022 年 10 月 19 日
OK, overlooked that in the indexing expression, yes.
Use
F(ix)=log(cons(ix));
F=F-l.^aux/aux;
instead.
I've some other commitments have to get ready for; if get a chance I'll see if I can find enough time to figure out what it is this is actually doing.
First guess that might help would be to look at meshgrid examples if what your doing is trying to evaluate a function over a 2D grid as it looks like might be...for debugging/testing, I'd also suggest cutting your array sizes down to something that fits on one screen in the command line window -- it's MUCH easier to figure out what's going on when the sizes are like 5x3 or somesuch instead of "cast of thousands" and the logic is no different with size.
採用された回答
dpb
2022 年 10 月 19 日
function F = myfun(l,aprime,a_val,z_val)
aux = 3.5;
cons = a_val+l*z_val - aprime;
F=-inf(size(cons));
indx=cons>0; % feasible!
F(indx)=log(cons(indx));
F=F-l.^aux/aux;
end
Using the above modfication of your function and a slightly smaller array size just so it wouldn't take too long I got the following results on local machine --
>> amm % original
Elapsed time is 2.460195 seconds.
>> amm % above modification
Elapsed time is 1.465042 seconds.
>>
That's about a 40% reduction which is not quite a factor of 2X, but getting close....
I didn't put the original back to test, but when I bumped the array sizes back to match yours, it was
>> amm
Elapsed time is 3.748369 seconds.
>>
Using the 60% factor, that would scale to about 6.3 sec for the original that is in ballpark of the above...
I did look at the code some more; it doesn't appear to be particularly amenable to further vectorizing.
I would suspect although I did not investigate further that the bigger speed up could be gained by screening the input vectors for the ranges that can't be feasible first and only iterating over those instead of the full, dead-ahead, just try 'em all approach.
その他の回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Resizing and Reshaping Matrices についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!