Is the pile-up of ones into the elements of a matrix parallelizable by parfor?

Hello,

Is there any way to do the following pile-up in matrix A using parfor?

A = zeros(1,J);
for i=1:n
  j = f(i); % f returns index j = 1,2,...J
  A(j) = A(j)+1;
end

Thanks in advance, Abi

[Jan]

Is "f" a function or a vector?

[Abi mehranian]

It's a time-consuming function!

Answer 1

[Matt J]

You don't really want the A(j)=A(j)+1 inside the PARFOR, because that's not a data parallel operation. Presumably f(i) is the hard computation, in which case you could do,

    A=zeros(J,1);
    j=ones(n,2);

    parfor i=1:n

            j(i,1) = f(i); % f returns index j = 1,2,...J

    end

    A=accumarray(j,1,size(A)).'

[Walter Roberson]

Right. And if f is simply a matrix of indices with no parallel computation in itself, then use

A = accumarray(f, 1, [1 J])

[Jan]

I do not agree. Imagine that f(1) replies [1,2] and f(2) replies [2,3]. Then A(j) = A(j) + 1 is correct and cannot be replaced by j(i,1) = f(i) and a following accumarray().

[Walter Roberson]

I am not sure what you mean by "reply" here?

You seem to be listing 2 values, but this is not a cell array.

[Jan]

@Walter: I meant that the function f() "replies" vectors of eventually different lengths. The OP used the term "returns" for this:

j = f(i); % f returns index j = 1,2,...J

What is not a cell array?

[Abi mehranian]

Many thanks for your comments,

As i mentioned f returns the index j, which takes values between 1,..J. the function f is the most time consuming part. As you suggested, A should be pile up out side of the parfor loop. actually, i was doing that by another for loop, which in my case should go through 200 millions loops. Thanks to Matt, i guess accumarray is a prudent alternative to that for loop. Here is a sample test (with f as a vector). Merci,

J = 20; n = 8000;

A=zeros(J,1);

j=ones(n,2);

f=floor(1+rand(1,n)*J);

parfor i=1:n

    j(i,1) = f(i); % f returns index j = 1,2,...J

end

A=accumarray(j,1,size(A)).'

[Jan]

So f replies a scalar? Then my suggestion to use a cell is not useful.

Answer 2

[Jan]

Assuming that f() is the most time-consuming part:

C = cell(1, n);
parfor iC = 1:n
  C{iC} = f(iC);
end

A = zeros(1,J);
for iC = 1:n
  index    = C{iC};
  A(index) = A(index) + 1;
end

Unfortunately I cannot test this, but collecting the results of f() in a cell should work, while A(j)=A(j)+1 cannot be performed inside a PARFOR due to a concurrent data access to the values of A.

[Abi mehranian]

Thanks Jan, you're right A can not be sliced between workers. as Matt and you suggested the parfor should be followed by accumarray in my case or a for loop in general.

[Matt J]

You could actually split the accumarray computation using PARFOR, but that means broadcasting multiple copies of A from the labs. For large J, I don't know if it would be worth it. My example below involves the cell-splitting utility MAT2TILES, available here

    j=ones(n,2);

    parfor i=1:n

            j(i,1) = f(i); % f returns index j = 1,2,...J

    end

    m=matlabpool('size');
    chunksize=ceil(J/m);
    jC = mat2tiles(j,[chunksize, inf]);

    parfor i=1:length(jC) 
     A=A+accumarray(jC{i},1,[J,1]);
    end

[Abi mehranian]

Thank you so much, that's fine too, but my 32 GB machine gets out of memory for J = 1e8, in any case, accumarray is much faster than a for-loop.