Vectorize the following loop

Hi all,

I'm trying to vectorize the following loop to speed it up

c = cumsum(weights);
A = ones(1,n);
x = rand(1,n);
for i = 1:n
     j = find(c > x(i) ,1,'first');
     A(i) = j;
end

where weights is an array of doubles which sum to 1 and n <= size(weights).

Any help would be grand!

B

Answer 1

[Oleg Komarov]

n = 3e4;
weights = rand(n,1);
c = cumsum(weights/sum(weights)); 
A = zeros(1,n); 
x = rand(1,n); 

tic
for i = 1:n 
    j = find(c > x(i) ,1,'first'); 
    A(i) = j; 
end
toc

% WARNING: only if weights are monotonically increasing

tic
[B1,B1] = histc(x,c);
toc

tic
B2 = n-sum(bsxfun(@gt,c,x))+1; % Goes fast in out of memory
toc

isequal(A,B1+1,B2) % Don't forget to add 1 to histc result

LOOP  : Elapsed time is 2.247998 seconds.
HISTC : Elapsed time is 0.005381 seconds.
BSXFUN: Elapsed time is 1.770875 seconds.

[Sean de Wolski]

Nice!

[Andrew Newell]

In recent versions of MATLAB you can use

[~,B1] = histc(x,c)

(but it doesn't make it any faster).

[Matt Fig]

Well-played, Oleg!

[Jan]

accepted by JSimon

Answer 2

[Sean de Wolski]

Note sure if it'll be faster but:

[row col] = find(bsxfun(@gt,c(:)',x(:)));
A2 = accumarray(row,col,[],@min)';

[Matt Fig]

Your intuition is correct. On my machine this is 50 times slower than the simple loop, using:

A = magic(1000);

weights = A(1,:)/sum(A(1,:));

n = 1000;

[Sean de Wolski]

It would only get slower as the matrices get bigger. I think Ben's elementary, but properly constructed, FOR-loop is probably optimal.

I just realized and am kind of surprised FIND doesn't have a dimensional argument.