vectorization of for loop

Hi,

I made a 2d matrix with two for loops:

for k = 1:32
for l = 1:32
P_new(l,k) = P_old(l) + (LODF(l,k) * P_old(k));
end
end

P_old is here a 32 x 1 matrix and LODF is a 32 x 32 matrix which is already computed. How can I vectorize this code to avoid the for loops? Thanks in advance.

Answer 1

[Matt J]

 P_new= bsxfun(@times, LODF, P_old.');
 P_new= bsxfun(@plus, P_new,P_old);

[Jip]

Although this code avoids the for loops, it is not faster, which is my purpose. Actually the code is much slower. Any other suggestions?

[Matt J]

For 32x32 data, I wouldn't be surprised if the for-loop was the fastest approach. For larger sizes, however, the vectorized approach will start to show superior performance, e.g.,

    N=3200;
    LODF=rand(N); P_old=rand(N,1);

    tic;
    P_new=zeros(size(LODF));
    for k = 1:N
    for l = 1:N
    P_new(l,k) = P_old(l) + (LODF(l,k) * P_old(k));
    end
    end
    toc;
    %Elapsed time is 0.102201 seconds.

    tic;
     P_new= bsxfun(@times, LODF, P_old.');
     P_new= bsxfun(@plus, P_new,P_old);
    toc
    %Elapsed time is 0.043591 seconds.

If you're not happy with the speed of your code, you should show us the slow part in its entirety. The small part you've shown is pretty fast, in and of itself.

Answer 2

[Jos (10584)]

for loops are pretty fast when you use pre-allocation

P_new = zeros(32,32) ; 
for k = 1:32
  for l = 1:32
    P_new(l,k) = P_old(l) + (LODF(l,k) * P_old(k));
  end
end

Answer 3

[Jip]

Yes I know, I already did that.