vectorizing a nested loop where one loop variable depends on the other

Question

I've recently learned how to vectorize a "simple" nested loop in a previous question I've asked. However, now I'm trying also to vectorize the following loop

A=rand(80,80,10,6,8,8);

I=rand(size(A1,3),1);
C=rand(size(A1,4),1);
B=rand(size(A1,5),1);

for i=1:numel(I)
    for v=1:numel(C)
        for j=1:numel(B)
            for k=1:j
                A(:,:,i,v,j,k)= A(:,:,i,v,j,k)*I(i)*C(v)*B(j)*((k-1>0)+1);               
            end
        end
    end
end

So now k depends in j... What have I tried so far: The combination of j and k terms (i.e. B(j)*((k-1>0)+1) gives a triangular matrix that I manage to vectorize independently:

  B2=tril([ones(8,1)*B']');
  B2(2:end,2:end)=2*B2(2:end,2:end);

But that gives me the (j,k) matrix properly and not a way to use it to vectorize the remaining loop. Maybe I'm in the wrong path too... So how can I vectorize that type of loop?

Community · Accepted Answer · 2017-05-23 12:13:54Z

13

In one of your comments to the accepted solution of the previous question, you mentioned that successive bsxfun(@times,..,permute..) based codes were faster. If that's the case, you can use a similar approach here as well. Here's the code that uses such a pattern alongwith tril -

B1 = tril(bsxfun(@times,B,[1 ones(1,numel(B)-1).*2]));
v1 = bsxfun(@times,B1, permute(C,[3 2 1]));
v2 = bsxfun(@times,v1, permute(I,[4 3 2 1]));
A = bsxfun(@times,A, permute(v2,[5 6 4 3 1 2]));

edited May 23, 2017 at 12:13

CommunityBot

11 silver badge

answered Oct 3, 2014 at 23:00

Divakar

222k19 gold badges273 silver badges374 bronze badges

Sign up to request clarification or add additional context in comments.

3 Comments

Max Over a year ago

wonderful! it's much more elegant and runs 25% faster than @natan's solution.

Divakar Over a year ago

@Max Awesome! Good to know that too!

rayryeng Over a year ago

This solution reminds me of Ramanujan. I have absolutely no idea how in hell you came up with that answer.

bla · Accepted Answer · 2014-10-03 07:08:07Z

You were close. The vectorization you proposed indeed follows the (j,k) logic, but doing tril adds zeros in places where the loop doesn't go into. Using the solution for your previous question (@david's) is not complete as it multiplies all elements including these zero value elements that the loop doesn't go into. My solution to that, is to find these zero elements and replace them with 1 (so simple):

Starting with your code:

B2=tril([ones(8,1)*B']');
B2(2:end,2:end)=2*B2(2:end,2:end);

and following the vectorization shown in the previous question:

s=size(A);
[b,c,d]=ndgrid(I,C,B2);
F=b.*c.*d;
F(F==0)=1; % this is the step that is important for your case.
A=reshape(A,s(1),s(2),[]);
A=bsxfun(@times,A,permute(F(:),[3 2 1]));
A=reshape(A,s);

For the size of A used in the question this cuts about 50% of the run time, not bad...

Collectives™ on Stack Overflow

vectorizing a nested loop where one loop variable depends on the other

2 Answers 2

3 Comments

Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

3 Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related