Problems in Parfor Implemenation

Hi,

I am new to parallel computing and I am trying to parallelize my (relatively intensive) code so as to make it run faster. However, I facing a lot of issues that I am not able to adress. The code is the following:

% Initial guess for policy function
  c_pol  = max(cmin, r*ones(S,1) * b_grid); 
  c_poli = c_pol;                           % next iteration

  % Guess for V_t+1 
  load CRRA.mat    % loads c_pol and n_pol from CRRA case in original GL code (check line before the start of terminal SS part)
  v_pol= (1/(1-bet))*(c_pol.^(1-gam)/(1-gam))+ pssi*(((1-n_pol).^(1-eta))/(1-eta));                                      
  v_poli = v_pol;  % next iteration                               
  
dif = 1;
  while dif > tol_pol
     ui = Pr * (((-v_pol).^(-alpha)).*(c_pol.^(-gam))); 
     ui = ui(:, b_grid >= -phi);                                                                                  
     v_u = Pr *((-v_pol).^(1-alpha)); 
     v_un = v_u(:,b_grid>= -phi);     
     v_con = v_u(:,b_grid < -phi); 
      
    for s=1:S      %% here i need parallelization
        c = ((1+r) * bet * ((v_un(s,:)).^(alpha/(1-alpha)) .* ui(s, :))).^ (-1/gam);  
        n = max(0, 1 - fac(s)*c.^gameta);                          
        b = b_grid(b_grid >= -phi) / (1+r) + c - theta(s)*n - z(s);         
        v = c.^(1-gam).*((1/(1-gam))) + pssi * (((1-n).^(1-eta))/(1-eta)) - bet* ((v_un(s,:)).^(1/(1-alpha))); 
     
        if b(1) > -phi                                    
            c_c = linspace(cl(s), c(1), Ic); % ***                                                                                     
            n_c = max(0, 1 - fac(s)*c_c.^gameta);                               
            b_c = -phi/(1+r) + c_c - theta(s)*n_c - z(s);                                                  
            vl = cl(1,1).^(1-gam).*((1/(1-gam))) + pssi * (((1-n_c).^(1-eta))/(1-eta))- bet* ((v_un(1,1)).^(1/(1-alpha)));
            v_c = linspace(vl(1),v(1),Ic); 
            b = [b_c(1:Ic-1), b];        
            c   = [c_c(1:Ic-1), c];          
            v = [ v_c(1:Ic-1), v];            
      end

      c_poli(s, :) = interp1(b, c, b_grid, 'linear', 'extrap'); % true c_t+1        
      v_poli(s,:) = interp1(b, v, b_grid, 'linear', 'extrap'); % true v_t+1          
      v_poli(v_poli>0)=-1e-6;
    end

    % check convergence
    c_poli = max(c_poli, cmin); 
    v_poli = min(v_poli, cmin);  
    dif1 = max(max(abs(c_poli - c_pol)));
    dif2 = max(max(abs(v_poli-v_pol)));
    dif = max([dif1,dif2]);
   
    % update
    c_pol = c_poli;
    v_pol = v_poli;

  end

Trying to use ''parfor'' I get the following error ''Unable to classify the variable 'v_poli' in the body of the parfor-loop'', which I don't know how to fix.

Any help would be highly appreciated

Answer 1

[Walter Roberson]

Change

      v_poli(s,:) = interp1(b, v, b_grid, 'linear', 'extrap'); % true v_t+1          
      v_poli(v_poli>0)=-1e-6;

to

      v_polis = interp1(b, v, b_grid, 'linear', 'extrap'); % true v_t+1
      v_polis(v_polis>0)=-1e-6;
      v_poli(s,:) = v_polis;

Otherwise you are trying to test all of v_poli() including the rows being created by other loop iterations.

[Panagiotis Veneris]

Thanks very much. That worked beautifully. However the code become even slower. To be precise, the whole codes goes like this

for it = 1:maxit
  % Update (government) budget constraint
  tau = (pr(1)*nu + r/(1+r)*B) / (1 - pr(1)); % labor tax
  z   = [nu, -tau*ones(1,S-1)];               % full transfer scheme (tau tilde in paper)

  % Find consumption at the lower bound of the state space: when both today
  % and tomorrow agents are at the borrowing limit ->constraint is binding
  % in both t and t+1
  fac = (pssi ./ theta) .^ (1/eta);
  for s = 1:S                        % loop over productivity states (13)
    cl(s) = fzero('find_cl',...      % objective
                  [cmin, 100],...    % interval for bisection
                  [],...             % options
                  s, -phi, -phi, r); % state, assets, interest rate
  end

......code attached in 1st message.

So, maybe I am not making efficient use of parallelization? Maybe I should use ''parfor'' somewhere else?

[Walter Roberson]

You could gain some performance improvement by storing the interp1() results into v_poli, but postponing the

      v_poli(v_poli>0)=-1e-6;

to after the parfor loop. I would not expect any major improvement though.

I recommend that you convert

          cl(s) = fzero('find_cl',...      % objective

to use a function handle to @find_cl . It looks like you are trying to pass in extra parameters, so convert that code according to http://www.mathworks.com/help/matlab/math/parameterizing-functions.html

[Panagiotis Veneris]

Could I gain performance by parallelizing the outer for loop ( ''for it = 1: maxit'')?

[Walter Roberson]

Unfortunately you have not happened to post enough of the code for me to assess whether that will work.

At the moment it is not clear to me which variables are updated when.

When you use parfor, the iterations must be independent of each other.

For example in the context of a budget cycle, the incoming debt might influence the budget process (because of interest payments), and it might hypothetically be the case that (for the purpose of discussion) that everything else is known variables except the incoming debt -- so (for discussion purposes) if you knew the incoming debt ahead of time for each year, perhaps you could process the years in parallel. But in practice the debt left over from one year is input to the next, so you would have to process the years in sequence instead of in parallel.

[Panagiotis Veneris]

% Housekeeping

%close all;

clear;

clc;

%=========================================================================

% 1. Set parameters

%=========================================================================

% Convenient to have parameters from labor supply equation as globals

global theta z fac gameta

% Fixed parameters

gam = 4; % inverse of EIS

frisch = 1; % avg Frisch elast.

r = 2.5/400; % ss interest rate

RA = 10; % risk aversion

alpha = 1-RA/gam; % a=0 (RA=gam) reduces to CRRA case

% Calibration targets

NE = 0.4; % avg hours of employed

nu_Y = 0.4; % UI benefit per GDP

B_4Y = 1.6; % liquid wealth per annual GDP

D1_4Y = 0.18; % HH debt to annual GDP in initial ss

D2_4Y = 0.08; % HH debt to annual GDP in terminal ss

% Leisure curvature follows right away

eta = 1/frisch * (1 - NE) / NE;

gameta = gam / eta;

% Numerical parameters

maxit = 500; % maximum number of iterations in calibration

cmin = 1e-6; % lower bound on consumption

tol_pol = 1e-10; % tolerance lvl for policy function

tol_dist = 1e-10; % tolerance lvl for distribution

tol_mkt = 1e-6; % tolerance lvl for market clearing

tol_cali = 1e-4; % tolerance lvl for calibration

% Asset grid

I = 200; % number of grid points for bonds

Ic = 100; % number of gridpoints for consumption % Ic =2;

bmin = -2; % lower bound

bmax = 50; % upper bound

b_grid = bmin + ((1:I)/I).^2 * (bmax - bmin); % denser for low values

db = 0.01; % step size for MPC

% Income shock process

load inc_process % loads x, Pr,pr

% 12-state Markov chain for employed states generated by Tauchen's method

% x : log productivity (real wage)

% Pr : transition matrix

% pr : invariant distribution

% Add unemployment

theta = [0; exp(x)];

S = length(theta); % number of states for productivity (13 states for productivity/real wage)

fin = 0.8820; % job-finding probability

sep = 0.0573; % separation probability

% new transition matrix

Pr = [1-fin, fin*pr; sep*ones(S-1, 1), (1-sep)*Pr];

% find new invariate distribution

pr = [0, pr];

dif = 1;

while dif > tol_dist

pri = pr*Pr;

dif = max(abs(pri-pr));

pr = pri;

end

% Initial guesses for calibrated parameters, NE ~ Y

bet = 0.8^(1/4); % discount factor

nu = nu_Y * NE; % UI benefits

B = B_4Y * NE * 4; % net supply of bonds

phi1 = D1_4Y * NE * 2; % borrowing constraint in initial ss = 0.1440

phi2 = D2_4Y * NE * 2; % borrowing constraint in terminal ss = 0.0640

pssi = NE^(-gam) * (1-NE)^eta; % disutility from labor as if representative agent

% Set rerun_initial = 1 for immediate calibration.

%rerun_initial = 1;

%rerun_terminal = 1;

%=========================================================================

% 2. Calibrate on Initial Steady State

%=========================================================================

%if rerun_initial == 1

% load par;

%end

% Set relevant borrowing constraint

phi = phi1;

% Preallocation

cl = ones(1, S); % consumption at borrowing constraint

n_pol = zeros(S, I); % labor policy

y_pol = zeros(S, I); % production policy

b_pol = zeros(S, I); % savings policy

mpcs = zeros(S, I); % MPCs

tic

disp('Computing initial equilibrium...');

parfor it = 1:maxit

% Update (government) budget constraint

tau = (pr(1)*nu + r/(1+r)*B) / (1 - pr(1)); % labor tax

z = [nu, -tau*ones(1,S-1)]; % full transfer scheme (tau tilde in paper)

% Find consumption at the lower bound of the state space: when both today

% and tomorrow agents are at the borrowing limit ->constraint is binding

% in both t and t+1

fac = (pssi ./ theta) .^ (1/eta);

for s = 1:S % loop over productivity states (13)

cl(s) = fzero('find_cl',... % objective

[cmin, 100],... % interval for bisection

[],... % options

s, -phi, -phi, r); % state, assets, interest rate

end

% A) Solve for consumption policy

%----------------------------------

% Initial guess for policy function

c_pol = max(cmin, r*ones(S,1) * b_grid); % guess for c_t+1

c_poli = c_pol; % next iteration

% Guess for V_t+1

load CRRA.mat % loads c_pol and n_pol from CRRA case in original GL code (check line before the start of terminal SS part)

v_pol= (1/(1-bet))*(c_pol.^(1-gam)/(1-gam))+ pssi*(((1-n_pol).^(1-eta))/(1-eta)); % guess for V_t+1 using c_pol, n_pol from CRRA case

%v_pol= (c_pol.^(1-gam)/(1-gam))+ pssi*(((1-n_pol).^(1-eta))/(1-eta)); % works fine as above

v_poli = v_pol; % next iteration

dif = 1;

while dif > tol_pol

% expected marginal utility tomorrow

ui = Pr * (((-v_pol).^(-alpha)).*(c_pol.^(-gam)));

ui = ui(:, b_grid >= -phi);

v_u = Pr *((-v_pol).^(1-alpha));

v_un = v_u(:,b_grid>= -phi);

v_con = v_u(:,b_grid < -phi);

%parfor s=1:S

for s=1:S

% unconstrained

c = ((1+r) * bet * ((v_un(s,:)).^(alpha/(1-alpha)) .* ui(s, :))).^ (-1/gam); % c_t

n = max(0, 1 - fac(s)*c.^gameta); % n_t

b = b_grid(b_grid >= -phi) / (1+r) + c - theta(s)*n - z(s); % b_t

v = c.^(1-gam).*((1/(1-gam))) + pssi * (((1-n).^(1-eta))/(1-eta)) - bet* ((v_un(s,:)).^(1/(1-alpha))); % v_t

% constrained

if b(1) > -phi

c_c = linspace(cl(s), c(1), Ic);

n_c = max(0, 1 - fac(s)*c_c.^gameta);

b_c = -phi/(1+r) + c_c - theta(s)*n_c - z(s);

vl = cl(1,1).^(1-gam).*((1/(1-gam))) + pssi * (((1-n_c).^(1-eta))/(1-eta))- bet* ((v_un(1,1)).^(1/(1-alpha))); % value function at the lower bound

v_c = linspace(vl(1),v(1),Ic); % value function (v_t) of constrained

b = [b_c(1:Ic-1), b]; % aggregate b_t

c = [c_c(1:Ic-1), c]; % aggregate c_t

v = [ v_c(1:Ic-1), v]; % aggregate v_t

end

c_poli(s, :) = interp1(b, c, b_grid, 'linear', 'extrap'); % true c_t+1 >0 everywhere

v_poli(s,:) = interp1(b, v, b_grid, 'linear', 'extrap'); % true v_t+1 <0 everywhere

v_poli(v_poli>0) = -1e-6;

% PARALLELIZED SECTION

%c_polis = interp1(b, c, b_grid, 'linear', 'extrap');

%c_poli(s,:) = c_polis;

%v_polis = interp1(b, v, b_grid, 'linear', 'extrap'); % true v_t+1 %

%v_polis(v_polis>0)=-1e-6;

%v_poli(s,:) = v_polis;

end

% check convergence

c_poli = max(c_poli, cmin); % enforce nonnegative constraint on consumption

v_poli = min(v_poli, cmin); %%%%CHECK!!! maybe we should have min

dif1 = max(max(abs(c_poli - c_pol)));

dif2 = max(max(abs(v_poli-v_pol)));

dif = max([dif1,dif2]);

% update

c_pol = c_poli;

v_pol = v_poli;

end

% Save other policy functions and MPCs

%parfor s = 1:S % parallelize

for s = 1:S

n_pol(s, :) = max(0, 1 - fac(s)*c_pol(s, :).^gameta);

y_pol(s, :) = theta(s) * n_pol(s, :);

b_pol(s, :) = max((1+r) * (b_grid + y_pol(s, :) - c_pol(s,:) + z(s)), -phi);

mpcs(s, :) = (interp1(b_grid, c_pol(s,:), b_grid + db, 'linear', 'extrap') - c_pol(s,:)) / db;

end

% B) Find invariant distribution (Iterate on prob. trans. matrix to find invariant distribution of bonds)

%--------------------------------------------------------------------------------------------------------

% Assign weights to adjacent grid points proportionally to distance

[~, ib_pol] = histc(b_pol, b_grid);

ib_pol(ib_pol<1) = 1;

wei = (b_pol - b_grid(ib_pol)) ./ (b_grid(ib_pol+1) - b_grid(ib_pol));

% Iterate asset transition matrix starting from uniform distribution

dif = 1;

pd = ones(S,I) / (S*I); % asset distribution (uniform) in period t, s: product. states, I: income states

while dif > tol_dist

pdi = zeros(S, I); % initialize period t+1 asset distribution

for s = 1:S % loop over productivity grid in period t

for i = 1:I % loop over asset grid in period t

for si = 1:S % loop over productivity grid in period t+1

pdi(si, ib_pol(s, i)) = (1 - wei(s, i)) * Pr(s, si) * pd(s, i) + pdi(si, ib_pol(s, i));

pdi(si, ib_pol(s, i) + 1) = wei(s, i) * Pr(s, si) * pd(s, i) + pdi(si, ib_pol(s, i) + 1);

end

end

end

% check convergence

dif = max(max(abs(pdi - pd)));

% make sure that distribution integrates (sums) to 1

pd = pdi / sum(sum(pdi));

end

% C) Check market clearing and calibration

%-------------------------------------------------

% Bond market clearing, i refers to current iteration

Bi = sum(sum(pd) .* b_grid);

res_mkt = abs(B - Bi);

% Calibration statistics, i refers to current iteration

Yi = sum(sum(pd .* y_pol)); % GDP

NEi = sum(sum(pd.*n_pol.*(n_pol>0))) / sum(sum(pd .* (n_pol>0))); % avg hours of employed

B_4Yi = Bi / Yi / 4; % debt ratio

Di = -sum(sum(pd) .* min(b_grid, 0));

D_4Yi = Di / Yi / 4; % debt ratio

nu_Yi = nu / Yi; % UI benefit ratio

res_cali = max(abs([B_4Yi, D_4Yi, nu_Yi, NEi] - [B_4Y, D1_4Y, nu_Y, NE]));

% Report convergence

disp(['Iteration ', num2str(it)]);

disp(['Bond mkt clearing: ', num2str(B - Bi)]);

disp(['Liquid wealth: ', num2str(B_4Yi - B_4Y)]);

disp(['Debt to GDP: ', num2str(D_4Yi - D1_4Y)]);

disp(['UI benefits: ', num2str(nu_Yi - nu_Y)]);

disp(['Avg hours: ', num2str(NEi - NE)]);

disp('-----------------------------------------------')

% Check convergence of both, update if necessary

if (res_cali < tol_cali) && (res_mkt < tol_mkt)

break;

else

% discount factor

btemp = -log(1-bet);

btemp = btemp - .1*(Bi - B_4Y*4*Yi);

bet = 1 - exp(-btemp);

% rest will be updated based on these

phi_d = phi * D1_4Y / D_4Yi;

nu_d = nu_Y * Yi;

pssi_d = pssi * ((1-NE) / (1-NEi))^eta;

phi = phi + 0.1 * (phi_d - phi);

nu = nu + 1 * (nu_d - nu);

pssi = pssi + 1 * (pssi_d - pssi);

% Update aggregates

B = B + 0.1 * (Bi - B); % net supply of bonds

end

end

% Aggregate statistics

Y1 = sum(sum(pd .* y_pol)); % GDP

C1 = sum(sum(pd .* c_pol)); % consumption

D1 = -sum(sum(pd) .* min(b_grid, 0)); % debt

N1 = sum(sum(pd .* n_pol)); % labor supply (counting unemployed)

MPC1 = sum(sum(pd .* mpcs)); % MPC

% Save results

phi1 = phi;

r1 = r;

pd1 = pd;

c_pol1 = c_pol;

b_pol1 = b_pol;

v_pol1 = v_pol;

% Figure 1 in paper

figure;

ii = (b_grid> -phi) & (b_grid < 50*Y1); % domain

subplot(1, 2, 1)

plot(b_grid(ii)/(4*Y1), c_pol(2, ii), b_grid(ii)/(4*Y1), c_pol(8, ii), '--', 'LineWidth', 1.3);

title('consumption');

box on; grid on;

% set(gca,'FontSize', 14);

axis([-2 13 0 .6])

subplot(1, 2, 2)

plot(b_grid(ii)/(4*Y1), n_pol(2, ii), b_grid(ii)/(4*Y1), n_pol(8, ii), '--', 'LineWidth', 1.3);

title('labor supply');

legend('\theta^2','\theta^8');

box on; grid on;

% set(gca,'FontSize', 14);

axis([-2 13 -.1 .7])

set(gcf,'Position', [440 378 700 300]);

set(gcf,'PaperPosition', [0 0 14 6]);

set(gcf,'PaperSize', [14 6]);

fig = gcf;

saveas(fig, 'fig1-policies.pdf');

This is the whole script. It would be quite useful to understand when (and where) i could gain from parallelization because it seems to me that the parts I have parallelized (see %parfor commands) lead to overheading, making the code almost 7 times slower (408 seconds vs 2661 seconds to execute).