Speed up vectorized code automatic expansion logical indexing
3 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Alessandro Maria Marco
am 19 Okt. 2022
Kommentiert: Alessandro Maria Marco
am 20 Okt. 2022
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 Kommentare
dpb
am 19 Okt. 2022
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.
Akzeptierte Antwort
dpb
am 19 Okt. 2022
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.
1 Kommentar
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Waveform Generation finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!