2014-10-02 174 views
6

我最近学会了如何向量化先前的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 ...什么都我想到目前为止: 的jk条款(即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)矩阵正确,而不是一种方法,用它来向量化剩余的循环也许我在错误的道路了。 ..那么我怎样才能矢量化这种类型的循环?

回答

11

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

太棒了!它更优雅,运行速度比@ natan的解决方案快25%。 – Max 2014-10-06 17:24:35

+0

@Max太棒了!很高兴知道这一点! – Divakar 2014-10-06 17:30:45

+6

这个解决方案让我想起[Ramanujan](https://en.wikipedia.org/wiki/Srinivasa_Ramanujan)。我完全不知道你是怎么想出这个答案的。 – rayryeng 2015-11-25 22:22:06

2

你很近。你提出的矢量化确实遵循(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%......