How to accumulate to a vector-indexed array?

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?

Answer 1

[Kelly Kearney]

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)

[David]

Kelly: Thanks, but it doesn't seem to work for me, at least in the original context. But I just realized, I'm still running R2007a, so possibly it's been fixed in a later version.

[Kelly Kearney]

No, Matlab from the very beginning allowed for indexing like this (certainly since v5, circa late 1990s).

Perhaps you could give a small example of your data, and the results you get from each of the above methods?

[David]

OK, here's a modified version of Kelly's example, a bit more like my case:

   clear

   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 vectorized method

   L2 = L;
   L2(K(1:njk)) = L2(K(1:njk)) + R(1:njk);

   isequal(L1, L2)

When I run this, L1 and L2 are not equal. And I can see in the workspace view that the max and min values are different.

Thanks in advance on any clues on why the difference.

[Kelly Kearney]

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)