Convert Nested for-Loops to parfor

I want to speed up my script and try to replace nested for loop with parfor. I know that I can parallelize either of the nested loops, but I cannot run both in parallel.

I am not sure how can I use a temporary arrays here if I need to save only recovery matrix.

Thanks for any help!

That is the part of my script I am working on:

%%
resultX0 = zeros(length(n), length(q), ntrials, M, N);
resultY0 = zeros(length(n), length(q), ntrials, M, N);
resultG = zeros(length(n), length(q), ntrials, M, N);
for i=1:length(n)
  for j=1:length(q)
      for k=1:ntrials
          [G,X0,Y0]=matrix_plantsubm(M,N,c,n(i),p,q(j));
          resultX0(i,j,k,:,:)=X0;
          resultY0(i,j,k,:,:)=Y0;
          resultG(i,j,k,:,:)=G;
      end
  end
end

%%
recovery = zeros(length(n), length(q));
error = zeros(length(n), length(q));
resultX=zeros(length(n),length(q), ntrials, M, N);
resultY=zeros(length(n),length(q), ntrials, M, N);

parfor i=1:length(n)
    for j=1:length(q)
        % Set gamma as in paper.
        gamma = 6/((q(j) - p)*n(i));
        
        for k=1:1:ntrials
            tic 
            [X,Y,Q, iter] = ADMMcopy(squeeze(resultG(i,j,k,:,:)),c, n(i),q(j), gamma,tau, opt_tol, verbose, maxiter);
            resultX(i,j,k, :, :)=X;
            resultY(i,j,k, :, :)=Y;
            C=resultX-resultX0;
a=arrayfun(@frobeniusnorm,C).^2;
b=arrayfun(@frobeniusnorm,resultX0).^2;
if a(i,j,k)/b(i,j,k)<opt_tol
                recovery(i,j)=recovery(i,j)+1;
                error(i,j)=error(i,j)+a(i,j,k)/b(i,j,k);
            else
                recovery(i,j)=0;
            toc
            save('Sparse_case.mat','recovery', '-v7.3');
        end
    end
    end
end

Answer 1

[Walter Roberson]

inside the parfor

trec = zeroes(1,length(q))

inside the for k

trec(1,j) = trec(1,j) + whatever

After the for j terminates

recovery(i,:)=trec;

You probably need something similar for error .

It would be better to rename error to avoid interfering with the important infrastructure function named error

Answer 2

[Matt J]

You can convert the first set of nested loops to a single parfor loop as below. The second set of nested-loops is non-parallelizable because the computation of C depends on all preceding iterations of the loops.

sz=[length(n), length(q)]; 
S=prod(sz);

resultX0 = zeros(S,ntrials, M, N);
resultY0 = zeros(S,ntrials, M, N);
resultG = zeros(S,ntrials, M, N);

parfor s=1:S
    
       [i,j] = ind2sub(sz,s);
    
          [G,X0,Y0]=matrix_plantsubm(M,N,c,n(i),p,q(j));
          resultX0(s,:,:,:)=repmat(X0,1,ntrials,1,1);
          resultY0(s,:,:,:)=repmat(Y0,1,ntrials,1,1);
          resultG(s,:,:,:)=repmat(G,1,ntrials,1,1);
 end
  
resultX0 = reshape(resultX,[sz ntrials M N]);
resultY0 = reshape(resultY0,[sz ntrials M N]);
resultG = reshape(resultG,[sz ntrials M N]);