Code vectorization problem

Hi

I have a piece of code with a for loop that I am trying to vectorize, but cannot find any easy way to do it. The code is

NG = 1e4;
N = 1e6;

a = floor(NG*rand(N, 1)) + 1;
b = rand(N, 1);
nn = sparse(NG, N);

for k = 1 : N

      i = a(k);
      nn(i, k) = b(k);

end

Does anyone have any suggestions? I tried to use linear indexing, but the matrix size is too large.

Thanks

Answer 1

[Paulo Silva]

nn=sparse(NG,N);
k=1:N;
nn(sub2ind([NG,N],a(k),k'))=b(k);

%Timing for NG = 10000 and N = 10000

%Elapsed time is 0.103679 seconds. -> with the for loop 
%Elapsed time is 0.003603 seconds. -> with the code of my answer

[Trevor]

Well it works, but unfortunately it is about 20-40 times slower than the current for-loop version

[Paulo Silva]

try it now, creating sparse arrays seems to be faster than arrays full of zeros

[Trevor]

Ok, now I time the vectorized code as being faster. On my computer I get:

(1) 0.068 sec with the for-loop

(2) 0.048 sec with the vectorized code

The speedup though is not quite as high, but it works, thanks. I'll wait for a day or so to see if anyone else has any ideas, otherwise I will accept this answer. I was hoping to get some ideas from this problem to speedup another for-loop, but I don't think the same trick will work. The other for-loop is

nn = zeros(NG, 1);

for k = 1 : N

i = a(k);

nn(i) = nn(i) + b(k);

end

I had thought that if I could form an NGxN matrix and vectorize that (like your answer has shown), then by summing the columns I could have a vectorized form of the new for-loop above. However the new for-loop above is much faster than the vectorized NGxN code. Any ideas how I might vectorize the new for-loop?

[Paulo Silva]

nn=accumarray(a, b); %no need to pre-allocate nn or create k

Answer 2

[Teja Muppirala]

You never want to build a sparse using a loop like this.

The SPARSE command is designed to handle that entire loop straight from the command line:

nn = sparse(a,1:N,b,NG,N);

Comparing this with a loop, you can see that calling the SPARSE command with the correct arguments works orders of magnitude faster.

%%%%6 seconds
NG = 1e4;
N = 1e5;
a = floor(NG*rand(N, 1)) + 1;
b = rand(N, 1);

tic
nn = sparse(NG, N);
for k = 1 : N
      i = a(k);
      nn(i, k) = b(k);
end
toc

%%%%0.007 seconds
tic
nn2 = sparse(a,1:N,b,NG,N);
toc

%%%%They give equal answers
isequal(nn,nn2)

[Trevor]

I originally used the SPARSE command in this way, but as I said in the comment to the other answer, by writing my question in this way I was hoping to get ideas to solve the second for-loop mentioned in that comment. SPARSE also works on that for-loop, but is much slower. Anyway, ACCUMARRAY works and is almost as fast as the second for-loop. Thanks for replying though.