Program runs normally in Matlab but won't work in Parallel Toolbox

Question

I recently asked about a related problem here:

https://stackoverflow.com/questions/21171836/storing-matrix-output-in-higher-dimensional-matrix

I now simply want to run the code in a parallel configuration, but when I do I get the following error:

??? Error using ==> parallel_function at 598 Error in ==> OIRE at 136 Undefined function or variable "best_index".

Error in ==> OIRE_MSE_test at 73 parfor t=1:nsims

Error in ==> OIRE_MSE_test at 95 [b_OIRE OIRE_opt_b(:,:,t)]=OIRE(y,x,iter);

The code works fine so long as I drop the matlab pool open/close and the parfor command. How come this won't run in parallel toolbox?

clc;
n=100;
p=7; 
alpha1=0.999999;
error_vol=0.1;

iter=1000;
nsims=200;

OIRE_opt_b=zeros(7,3,nsims);
OIRE_opt_MSE=zeros(1,3,nsims);
OIRE_opt_index=zeros(1,3,nsims);

GIREI_opt_b=zeros(7,3,nsims);
GIREI_opt_MSE=zeros(1,3,nsims);
GIREI_opt_index=zeros(1,3,nsims);

GIREII_opt_b=zeros(7,3,nsims);
GIREII_opt_MSE=zeros(1,3,nsims);
GIREII_opt_index=zeros(1,3,nsims);

LRRE_opt_b=zeros(7,3,nsims);
LRRE_opt_MSE=zeros(1,3,nsims);
LRRE_opt_index=zeros(1,3,nsims);



 matlabpool open


x=zeros(n,p);

for i=1:n
    z_i4=normrnd(0,1);
    x(i,p)=z_i4;
    for j=1:p-1
        x(i,j)=x(i,j)+alpha1*z_i4;
        x(i,j)= x(i,j)+(1-alpha1^2)^(0.5)*normrnd(0,1);
    end  
end

b_act=[5;1;10;-20;200;30;-2];

parfor t=1:nsims

    residuals=normrnd(0,error_vol,n,1);
    y=x*b_act + residuals;
    y_store(:,t)=y;
    y=y_store(:,t);

    [b_OIRE OIRE_opt_b(:,:,t) OIRE_opt_MSE(:,:,t) OIRE_opt_index(:,:,t)]=OIRE(y,x,iter);

end

called function

function [b_OIRE OIRE_opt_b OIRE_opt_MSE OIRE_opt_index]=OIRE(y,x,iter)  
dim=1;
pool=[10,100,1000,10000,1000000,10000000];
count=0;

[n, p]=size(x);

b=x\y; 
b_OIRE = b; % [#1] initialize b_LRRE as b
sigma_sq=((y-x*b)'*(y-x*b))/(n-p); %'

b_act=[1;0;1;1;0;1;1];

econFlag=0;
[U,sigma,V] = svd(x,econFlag);


U1=U(:,1:p);
d=zeros(p,1);
d=diag(d);

alpha=V'*b_OIRE; %'
Delta=sigma.^1;        
Delta=diag(Delta);
f=Delta.*alpha;
F=diag(f);
Theta=sum(f);
c=p^2*sigma_sq+p*Theta^2;
g=Theta*sum(alpha);
I=ones(p,1);
a=sigma_sq*I+Theta*f;
b=F*alpha;
k=zeros(p,1);
A=sigma_sq*eye(p)+F.^2;
varRho=(g-a'*pinv(A)*b)*pinv(c-a'*pinv(A)*a); 
k=pinv(A)*b-varRho*pinv(A)*a;
K=diag(k);
D=varRho*I*I';

b_OIRE= V*(K+D)*U1'*y;


MSE=(k'*A*k)+(2*varRho*a'*k)-(2*b'*k)+(c*varRho^2)-(2*g*varRho)+(alpha'*alpha);


best_OIRE_MSE=MSE;
best_b_OIRE=b_OIRE;


for jj=1:iter           % [## "iter" denotes the iteration number]

    alpha=V'*b_OIRE;        %'
    Delta=sigma.^1;         % [Error! not sigma.^2 but sigma.^1]
    Delta=diag(Delta);
    f=Delta.*alpha;
    F=diag(f);
    Theta=sum(f);
    c=p^2*sigma_sq+p*Theta^2;
    g=Theta*sum(alpha);
    I=ones(p,1);
    a=sigma_sq*I+Theta*f;
    b=F*alpha;
    k=zeros(p,1);
    A=sigma_sq*eye(p)+F.^2;
    varRho=(g-a'*pinv(A)*b)*pinv(c-a'*pinv(A)*a);
    k=pinv(A)*b-varRho*pinv(A)*a;
    K=diag(k);
    D=varRho*I*I';

    b_OIRE= V*(K+D)*U1'*y;

    MSE=(k'*A*k)+(2*varRho*a'*k)-(2*b'*k)+(c*varRho^2)-(2*g*varRho)+(alpha'*alpha);

    if(MSE<best_OIRE_MSE)

        best_b_OIRE=b_OIRE;
        best_OIRE_MSE=MSE;
        best_index=jj+1;
    end

    if( any(jj == pool))
        count=count+1;
        OIRE_opt_b(:,count)=best_b_OIRE;
        OIRE_opt_MSE(count)=best_OIRE_MSE;
        OIRE_opt_index(count)=best_index;
    end


end


end

matlabpool close

Answer 1

Your problem is that best_index is never defined anywhere other than in the MSE<best_OIRE_MSE if-statement.

This is what I suspect is happening. When you run your code without the parfor, it iterates through sequentially (i.e. jj==1 will always occur before jj==2 etc). Without looking through the intricacies of your code, I suspect this means that MSE<best_OIRE_MSE will always be true before (or on the same iteration as) any(jj==pool) is true. This means that by the time your code reaches the any(jj == pool) if-statement, best_index will always have been set.

The problem comes with the fact that parfor offers no guarantees as to the execution order of your loop. jj==10 might well run before jj==2. This means that there is no guarantee that best_index will be defined by the time it reaches the any(jj == pool) if-statement - hence your error.

My suggestion would be to define best_index outside your main parfor loop to something like -1, and then ignore any cases where OIRE_opt_index is negative.