How to accumulate to a vector-indexed array?
1 view (last 30 days)
Show older comments
I want to vectorize something in the following form:
for jk = 1:njk
k = K(jk);
L(k) = L(k) + R(jk);
end
where the range of the k values is less than that of the jk values, so some k values occur for multiple values of jk. The obvious "solution" would be
L(K(1:njk)) = L(K(1:njk)) + R(1:njk);
but this gives wrong results, either because (1) L(K(1:njk)) occurs on both the LH and RH sides, or because (2) vectored indexing cannot be used on the LH side (I'm not sure whether (1) or (2) is the reason). A possible approach might be if there is something equivalent to the "+=" operator in C, which would allow
L(K(1:njk)) += R(1:njk);
Also the vectors happen to be large (>10000 elements) so anything that involves creating an intermediate matrix ix not likely to be practical.
Is there a way to vectorize this?
0 Comments
Accepted Answer
Kelly Kearney
on 4 Apr 2014
I don't see anything wrong with the syntax you proposed...
L = rand(100,1);
R = rand(100,1);
K = randperm(100);
njk = 5;
% The loop method
L1 = L;
for jk = 1:njk
k = K(jk);
L1(k) = L1(k) + R(jk);
end
% The vectorized method
L2 = L;
L2(K(1:njk)) = L2(K(1:njk)) + R(1:njk);
% Yup, they're the same...
isequal(L1, L2)
4 Comments
Kelly Kearney
on 7 Apr 2014
Ah, I see... the differences are due to the repeated indices. When you assign using repeated indices in a single call, only the last instance is used. For example:
>> a = rand(5,1); a([1 2 1]) = [1 2 3]
a =
3
2
0.41407
0.89632
0.51201
I'm sure that must be documented somewhere, though I can't find it in my quick search.
In your case, I think accumarray is the best option. It's a slightly confusing function when you're first introduced to it, but exactly what's needed in this situation.
njk = 1000;
L = rand(1,100);
R = rand(1,1000);
K = int16(sort(rand(1,1000))*100 + 0.5); % values 1 to 100
% The loop method
L1 = L;
for jk = 1:njk
k = K(jk);
L1(k) = L1(k) + R(jk);
end
% The vecotorized method
L2 = L;
nk = max(K);
rsum = accumarray(K', R', [nk 1], @sum);
L2(1:nk) = L2(1:nk) + rsum';
% Are they the same? Well, close (some rounding error)...
isequal(L1, L2)
max(L1 - L2)
More Answers (0)
See Also
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!