vectorizing a nested loop where one loop variable

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?

Answer 1

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...

Answer 2

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]));