向量化嵌套循环，其中一个循环变量依赖于另一个

Question

我最近在我之前问过的question 中学习了如何对“简单”嵌套循环进行矢量化。但是，现在我也在尝试对以下循环进行矢量化

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

所以现在k 依赖于j...到目前为止我尝试了什么： j 和 k 术语的组合（即B(j)*((k-1>0)+1) 给出了一个我设法独立矢量化的三角矩阵：

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

但这给了我正确的 (j,k) 矩阵，而不是使用它来向量化剩余循环的方法。也许我也走错了路……那么我怎样才能矢量化这种类型的循环呢？

Answer 1

你很亲密。您提出的矢量化确实遵循 (j,k) 逻辑，但是执行 tril 在循环不进入的地方添加零。使用您之前的问题（@david's）的解决方案并不完整，因为它将所有元素相乘，包括循环未进入的这些零值元素。我的解决方案是找到这些零元素并将它们替换为 1（非常简单）：

从您的代码开始：

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

并遵循上一个问题中显示的矢量化：

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

对于问题中使用的A 的大小，这减少了大约 50% 的运行时间，还不错......

Answer 2

在one of your comments 对上一个问题的公认解决方案中，您提到基于bsxfun(@times,..,permute..) 的连续代码更快。如果是这种情况，您也可以在此处使用类似的方法。这是使用这种模式的代码以及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]));