Parallelize nested loops with parfor

5 visualizaciones (últimos 30 días)
Alessandro Maria Marco el 27 de En. de 2023
Comentada: Alessandro Maria Marco el 27 de En. de 2023
I'm trying to speed up the simulation of some panel data. I have to simulate first over individuals (for ii from 1 to N) and then for each individual over age from 1 to JJ. The code is slow because inside the two loops there is a bilinear interpolation to do. Since the iterations in the outer loop are independent, I tried to use parfor in the outer loop, but I get the error message "the parfor cannot run due to the way the variable hsim is used".
Could someone explain why and how to solve the problem if possible? Any help is greatly appreciated!
%% Generate fake data for MWE
clear;clc
Nsim = 500;
JJ = 66;
na = 50;
nh = 30;
nz = 7;
a_grid = linspace(-15,100,na)'; % NOT necessarily equispaced!
h_grid = linspace(0,20,nh)'; % NOT necessarily equispaced!
apol = rand(na,nh,nz,JJ);
hpol = rand(na,nh,nz,JJ);
z_sim_ind = unidrnd(nz,Nsim,JJ);
%---------- START REAL CODE
tic
% a_grid and h_grid are column vectors with 50 and 35 elements,
% respectively
% z_sim_ind is a matrix of integers with size [Nsim,JJ]
a_sim = zeros(Nsim,JJ);
h_sim = zeros(Nsim,JJ);
% Find point on a_grid corresponding to zero assets
aa0 = 8;%find_loc(a_grid,0.0);
% Zero housing
hh0 = 1;
a_sim(:,1) = a_grid(aa0);
h_sim(:,1) = h_grid(hh0);
for ii=1:Nsim %with for loop it works but is very slow. parfor is illegal
for jj=1:JJ-1
z_c = z_sim_ind(ii,jj);
apol_interp = griddedInterpolant({a_grid,h_grid},apol(:,:,z_c,jj));
hpol_interp = griddedInterpolant({a_grid,h_grid},hpol(:,:,z_c,jj));
a_sim(ii,jj+1) = apol_interp(a_sim(ii,jj),h_sim(ii,jj));
h_sim(ii,jj+1) = hpol_interp(a_sim(ii,jj),h_sim(ii,jj));
end
end
toc
Elapsed time is 1.200044 seconds.
0 comentariosMostrar -2 comentarios más antiguosOcultar -2 comentarios más antiguos

Iniciar sesión para comentar.

Respuestas (2)

Matt J el 27 de En. de 2023
Editada: Matt J el 27 de En. de 2023
I don't think it makes sense to use any loops here at all. Just make vectorized calls to your griddedInterpolant objects.
[~,~,m,n]=size(apol);
[~,Jcol]=ndgrid(1:Nsim,1:JJ);
F=griddedInterpolant({a_grid,h_grid,1:m,1:n},apol);
a_sim=F(a_sim(:),h_sim(:), z_sim_ind(:),Jcol(:));
F.Values=hpol;
h_sim=F(a_sim(:),h_sim(:), z_sim_ind(:),Jcol(:));
a_sim=circshift( reshape(a_sim, Nsim,JJ) ,[0,1]);
h_sim=circshift( reshape(h_sim, Nsim,JJ ,[0,1]);
a_sim(:,1) = a_grid(aa0);
h_sim(:,1) = h_grid(hh0);
2 comentariosMostrar NingunoOcultar Ninguno
Alessandro Maria Marco el 27 de En. de 2023
Editada: Alessandro Maria Marco el 27 de En. de 2023
Matt J el 27 de En. de 2023
I think I fixed the discrepancy, but I would really need your input data (in a .mat file) to know for sure.

Iniciar sesión para comentar.

Matt J el 27 de En. de 2023
Editada: Matt J el 27 de En. de 2023
Never mind my other answer. I didn't notice that a_sim and h_sim were recursively defined. You can't avoid a loop, but you only need the loop over jj. It cannot be parallelized.
Nsim = 50000;
JJ = 66;
% a_grid and h_grid are column vectors with 50 and 35 elements,
% respectively
% z_sim_ind is a matrix of integers with size [Nsim,JJ]
a_sim = zeros(Nsim,JJ);
h_sim = zeros(Nsim,JJ);
% Find point on a_grid corresponding to zero assets
aa0 = 8;%find_loc(a_grid,0.0);
% Zero housing
hh0 = 1;
a_sim(:,1) = a_grid(aa0);
h_sim(:,1) = h_grid(hh0);
%%%%%%%%%%%%%%%%%%%%%%%
k=(1:Nsim)';
apol_interp = griddedInterpolant({a_grid,h_grid,k},apol(:,:,k,1));
hpol_interp = griddedInterpolant({a_grid,h_grid,k},hpol(:,:,k,1));
for jj=1:JJ-1
z_c = z_sim_ind(:,jj);
I=sub2ind( [Nsim,JJ], zc, repelem(jj,Nsim,1) );
apol_interp.Values = apol(:,:,I);
hpol_interp.Values = hpol(:,:,I);
a_sim(:,jj+1) = apol_interp(a_sim(:,jj),h_sim(:,jj),k);
h_sim(:,jj+1) = hpol_interp(a_sim(:,jj),h_sim(:,jj),k);
end
2 comentariosMostrar NingunoOcultar Ninguno
Alessandro Maria Marco el 27 de En. de 2023
Thanks again but I get an error unfortunately. I edited my original question with a minimum working example, so users can run the code themselves and check the execution time.
Alessandro Maria Marco el 27 de En. de 2023
The problem is that apol(:,:,k,1) gives an error because the 3rd dimension of the array "apol" has nz=7 elements, but k is a vector with Nsim>nz elements. I apologize if my original post was not clear enough. I added a MWE so everyhting should be clear now

Iniciar sesión para comentar.

Categorías

Más información sobre Data Type Identification en Help Center y File Exchange.

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by