Nested loop vectorization where one loop variable depends on another

Question

I recently learned how to vectorize a “simple” nested loop in a previous question that I asked. However, now I am trying to vectorize the next loop as well

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 on j ... What I have tried so far: The combination of the terms j and k (i.e. B(j)*((k-1>0)+1) gives a triangular matrix, which I manage to independently vectorize:

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

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

Answer 1

In one of your comments on the decision made on the previous question, you mentioned that the sequential bsxfun(@times,..,permute..) codes bsxfun(@times,..,permute..) were faster. In this case, you can also use a similar approach. Here's the code that uses such a template along with 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])); 

Answer 2

You were close. You the proposed vectorization does follow the logic (j, k), but tril adds zeros in places where the loop does not enter. Using the solution for your previous question (@ david's) is not complete, as it multiplies all elements, including these zero-value elements, which the loop does not include. My solution to this is to find these null elements and replace them with 1 (so simple):

Starting from your code:

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

and after 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 size A used in the question, this reduces about 50% of the execution time, not bad ...