Code vectorization problem

2011 7 5
2 个回答
回答已采纳
3 次查看（30 天）

0 个投票

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

请先登录，再回答此问题。

 采纳的回答

Paulo Silva
Paulo Silva 2011-7-5

0 个投票

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

4 个评论
显示 2更早的评论 隐藏 2更早的评论

Trevor
Trevor 2011-7-5
Well it works, but unfortunately it is about 20-40 times slower than the current for-loop version
Paulo Silva
Paulo Silva 2011-7-5
try it now, creating sparse arrays seems to be faster than arrays full of zeros
Trevor
Trevor 2011-7-5
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
Paulo Silva 2011-7-5
nn=accumarray(a, b); %no need to pre-allocate nn or create k

请先登录，再进行评论。

更多回答（1 个）

Teja Muppirala
Teja Muppirala 2011-7-5

0 个投票

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)

1 个评论
显示 -1更早的评论 隐藏 -1更早的评论

Trevor
Trevor 2011-7-5
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.

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by