MATLAB: How to modify this nested loop to work as a parfor without requiring a ridiculously large array

Question

I am trying to nest several for loops in a parfor loop. I have even read the documentation and several other queries/replies before posting here.

I have a large dataset that I need to iterate over, calculating a property which would result in a prohibitively large array if I sent each answer to a separate element, and doing this on one cpu will take a prohibitively long time.

Instead, what I am doing is reducing the array by another property, and then combining the calculated results in bins of this second property (essentially making a histogram).

The parfor section of my code looks like this:

 parfor i=1:box(1);
     for j=1:box(2);
         for k=1:box(3);
             for l=1:box(1);
                 for m=1:box(2);
                     for n=1:box(3);
 prop(horrendousfunction(i,j,k,l,m,n)) = prop(horrendousfunction(i,j,k,l,m,n)) + data(i,j,k)*data(l,m,n);
                     end
                 end
             end
         end
     end
 end

Trialling this on one cpu over i=j=k=l=m=n=[1:15] works fine and takes a few minutes.

The data array is the initial large array, but as written the code will iterate over every element of it many times within each parfor step, and therefore the data transmission I/O overhead shouldn't be too onerous with respect to the overall computation time. (Expected complete running time is ~1 month on 1 cpu, and I have several of these to do)

A few other (small, negligible I/O) arrays are also required.

horrendousfunction is essentially a way of mapping (i,j,k,l,m,n) into a single index h which has been previously split into bins earlier in the code. I will eventually want to plot(h,prop).

Now, after spending some time trawling the documentation and various previous questions on the topic, I realise that my initial efforts contravene the basic parfor assumption that j,k,l,m, and n (as components of the index of prop) should be constant within the parfor loop. I see that if I wanted to flesh out a behemoth array A(i,j,k,l,m,n) I could do so by making dummy variables and then combining it all with a command = A(i,:,:,:,:,:) at the end – but this will create the stupidly large array that I wish to avoid (I don't have enough HDD to store it).

Simply calculating trick = horrendousfunction(i,j,k,l,m,n) within the loop is also not an option, because I can't use trick as the index of prop either, since its value changes within the loop too.

Predetermining the complete mapping of [i,j,k,l,m,n] -> h is also not an option, since that would require a behemoth data array of its own, and require said array to be I/O'd to all the matlabpool workers. (And calculating it would take about as long as the original problem anyway) Also, the mapping is context-dependent, so I can't use the same mapping across the related family of problems I wish to solve.

Is there another way I can write the internals of this loop to take advantage of MATLAB's parallelism without requiring a behemoth output array?

Answer 1

I think I might use SPMD for this, instead of parfor, as below.

So the first modification to make is to abandon 6-dimensional indexing i,j,k,l,m,n. Replace subscripts (i,j,k) everywhere at the top level by a linear index u and (l,m,n) by a linear index v, so that your histogram update reduces simply to

prop(horrendousfunction(u,v)) = ...
    prop(horrendousfunction(u,v)) + data(u)*data(v);

You can call ind2sub() inside horrendousfunction as needed with marginal overhead, since it's already horrendous.

The next step is to write a modified and vectorized version of horrendousfunction(). Namely you will write

   h=horrendousVectorized(U,v)

to take a vector U of linear indices u and a scalar linear index v and return the appropriate vector of h indices, i.e., h(k)=horrendousfunction(U(k),v) for all k. This should be easy if, as you say,the calculation of h just involves abs() and simple arithmetic. All of those ops are nicely vectorizable for fixed v.

Then it's all just,

    B=prod(box);
    p=zeros(N,1);
    D=distributed(data(:));
    U=distributed((1:B).');
    for v=1:B
        datav=data(v);
        spmd
            h=horrendousVectorized(U,v);
            vals=D*datav;
            p=p+accumarray(h,vals,[N,1]);
        end
    end
   spmd
    p=gplus(p,1);
   end
   prop=gather(p{1});