Array indexing question for vectorization

I have three arrays A, B, and IND, with A containing some values that I want to update, B containing the values that I want to use to update A, and IND containing the indeces to update in A. IND may reference the same index in A several times. For example, I may have:

A = ones(1, 10);
B = 1:10;
IND = [1 1 2 5 5 5 5 6 8 9];

and the update I want to do is multiplication. I can achieve this by doing

for ii = 1:10
  A(IND(ii)) = A(IND(ii)) * B(ii);
end

which gives me my desired answer of

A =
   2     3     1     1   840     8     1     9    10     1

Is there a way I can do this in a vectorized operation, avoiding the for loop (I'm potentially doing this kind of operation thousands of times on large arrays)? Doing

A(IND) = A(IND) .* B

results in

A =
   2     3     1     1     7     8     1     9    10     1

Any tips greatly appreciated!

[Matt Fig]

Are A and B always like that? Could they be more like:

A = randi(N,1,M);
B = randi(K,1,J);

or not?

[Richard]

Hi Matt,

My initial problem specification was indeed lacking - see my reply to Jan below for a more accurate representation of the problem...

Answer 1

[Honglei Chen]

You can try the following trick

c = unique(IND)
d = accumarray(IND',B',[],@prod)
A(c) = A(c).*d(c)'

[Richard]

Fantastic, thanks! I didn't know about accumarray(), looks very useful.

[José-Luis]

This is not a good idea though. A for loop will be much faster in this case.

unique() is expensive

accumarray() is expensive

1. Loop: Elapsed time is 0.000443 seconds.
2. Accumarray: Elapsed time is 0.021888 seconds.

[Matt J]

There's no real need for UNIQUE, though, as far as I can see

 z=accumarray(IND(:),B,[10,1],@prod,1).'; 
 A=A.*z;

[Jan]

@Jose-Luis: Please post this as an answer such that I can vote for it. @Richard: Please post the typical size of the data. While the loop will have a nearly linear scaling, unique is less efficient dur to the internal sorting. If your data are sorted already, there are smarter versions than unique for the first part. accumarray suffers from calling Matlab from a Mex function.

[José-Luis]

True, but it still is slower than a for loop:

Elapsed time is 0.000003 seconds. %for

Elapsed time is 0.000828 seconds. %unique+accumarray

Elapsed time is 0.000515 seconds. %accumarray

[Richard]

Interesting, thanks. I'm still in the mindset of "for-loop bad", and aiming to replace them without necessarily thinking of the built-in function overheads...

Answer 2

[Jan]

Is your A pre-allocated? If IND is sorted, running the loop backwards wil be faster:

for ii = 10:-1:1
  A(IND(ii)) = A(IND(ii)) * B(ii);
end

If IND is not sorted, pre-allocate explicitly:

A = zeros(1, max(IND));

For more speed create a C-Mex function. If you are interested and have a compiler installed, I can post the few required lines of code.

[Richard]

Hi Jan,

Yes, A is preallocated. Actually the problem is a bit more complicated than I posted (note to self: make sure your example code is suitable!): in my actual application, A is a 2D matrix, and B is a 1D array with fewer elements than A, so the solution I've accepted doesn't quite work... e.g.

% A more accurate example...
A = zeros(20000, 10);
indR = randi(20000, 1, 500);         % row indeces
indC = randi(10, 1, 500);            % column indeces
IND = sub2ind(size(A1), indR, indC);
B = rand(1, 500);

tic
for ii = 1:length(IND)
  A(IND(ii)) = A(IND(ii)) * B(ii);
end
toc

(in this example with random indeces I think the instances of duplicate values in IND are quite low compared with the actual program). Thinking about it, my IND values could be sorted if I reorganised my data a bit prior to this (which may be worth the reorganisation overhead as this bit of the program is running in a simulation loop over many time steps)...

Thanks for the useful discussion everyone.

[Jan]

@Richard: A = zeros(...); ... A(k) = A(k) * ... leads to zeros. As in your original question your need ones().