如何向量化嵌套循环

Posted 2023-02-24

【中文标题】如何向量化嵌套循环

【英文标题】：How to Vectorize a Nested Loop

【发布时间】：2012-11-04 11:01:08

【问题描述】：

我无法想象如何矢量化这组循环。任何指导将不胜感激。

ind_1 = [1,2,3];
ind_2 = [1,2,4];
K = zeros(3,3,3,3,3,3,3,3,3);
pp = rand(4,4,4);

for s = 1:3
 for t = 1:3
  for k = 1:3
   for l = 1:3
    for m = 1:3
     for n = 1:3
      for o = 1:3
       for p = 1:3
        for r = 1:3
         % the following loops are singular valued except when
         % y=3 for ind_x(y) in this case
         for a_s = ind_1(s):ind_2(s)
          for a_t = ind_1(t):ind_2(t)
           for a_k = ind_1(k):ind_2(k)
            for a_l = ind_1(l):ind_2(l)
             for a_m = ind_1(m):ind_2(m)
              for a_n = ind_1(n):ind_2(n)
               for a_o = ind_1(o):ind_2(o)
                for a_p = ind_1(p):ind_2(p)
                 for a_r = ind_1(r):ind_2(r)
                  K(s,t,k,l,m,n,o,p,r) = K(s,t,k,l,m,n,o,p,r) + ...
                    pp(a_s, a_t, a_r) * pp(a_k, a_l, a_r) * ...
                    pp(a_n, a_m, a_s) * pp(a_o, a_n, a_t) * ...
                    pp(a_p, a_o, a_k) * pp(a_m, a_p, a_l);
                 end
                end
               end
              end
             end
            end
           end
          end
         end
        end
       end
      end
     end
    end
   end
  end
 end
end

编辑：

代码通过将每个索引的 pps 乘积的值相加一或两次来创建索引从 1 到 3 的 rank-9 张量，具体取决于 ind_1 和 ind_2 的值.

编辑：

这是一个 3d 示例，但请记住，pp 的索引只是简单置换的事实在 9d 版本中没有保留：

ind_1 = [1,2,3];
ind_2 = [1,2,4];
K = zeros(3,3,3);
pp = rand(4,4,4);

for s = 1:3
 for t = 1:3
  for k = 1:3
   % the following loops are singular valued except when
   % y=3 for ind_x(y) in this case
   for a_s = ind_1(s):ind_2(s)
    for a_t = ind_1(t):ind_2(t)
     for a_k = ind_1(k):ind_2(k)
      K(s,t,k) = K(s,t,k) + ...
        pp(a_s, a_t, a_r) * pp(a_t, a_s, a_k) * ...
        pp(a_k, a_t, a_s) * pp(a_k, a_s, a_t);
     end
    end
   end
  end
 end
end

【问题讨论】：

您能详细说明一下这段代码在做什么并添加 cmets 吗？

您能创建一个 2D 或 3D 示例来说明吗？与我们合作会更容易，创建示例可能会帮助您自己制定策略。

@tmpearce：已编辑以提供 3d 示例。

这是一些令人印象深刻的字母水果圈汤。

您是否查看过内置函数：mathworks.nl/help/matlab/ref/kron.html？如果可以使用此功能，我怀疑其他任何东西都会显示出更好的性能。

【参考方案1】：

哇！非常简单的解决方案，但并不容易找到。顺便说一句，我想知道你的公式是从哪里来的。

如果您不介意暂时丢失一点内存（2 次 4^9 数组 vs 3^9 之前），您可以在最后推迟第 3 和第 4 超平面的累积。

在 unix 机器上使用 octave 3.2.4 进行测试，它从 23s (67Mb) 下降到 0.17s (98Mb)。

function K = tensor9_opt(pp)

  ppp = repmat(pp, [1 1 1 4 4 4 4 4 4]) ;
  % The 3 first numbers are variable indices (eg 1 for a_s to 9 for a_r)
  % Other numbers must complete 1:9 indices in any order
  T = ipermute(ppp, [1 2 9 3 4 5 6 7 8]) .* ...
      ipermute(ppp, [3 4 9 1 2 5 6 7 8]) .* ...
      ipermute(ppp, [6 5 1 2 3 4 7 8 9]) .* ...
      ipermute(ppp, [7 6 2 1 3 4 5 8 9]) .* ...
      ipermute(ppp, [8 7 3 1 2 4 5 6 9]) .* ...
      ipermute(ppp, [5 8 4 1 2 3 6 7 9]) ;

  % I have not found how to manipulate 'multi-ranges' programmatically. 
  T1 = T (:,:,:,:,:,:,:,:,1:end-1) ; T1(:,:,:,:,:,:,:,:,end) += T (:,:,:,:,:,:,:,:,end) ;
  T  = T1(:,:,:,:,:,:,:,1:end-1,:) ; T (:,:,:,:,:,:,:,end,:) += T1(:,:,:,:,:,:,:,end,:) ;
  T1 = T (:,:,:,:,:,:,1:end-1,:,:) ; T1(:,:,:,:,:,:,end,:,:) += T (:,:,:,:,:,:,end,:,:) ;
  T  = T1(:,:,:,:,:,1:end-1,:,:,:) ; T (:,:,:,:,:,end,:,:,:) += T1(:,:,:,:,:,end,:,:,:) ;
  T1 = T (:,:,:,:,1:end-1,:,:,:,:) ; T1(:,:,:,:,end,:,:,:,:) += T (:,:,:,:,end,:,:,:,:) ;
  T  = T1(:,:,:,1:end-1,:,:,:,:,:) ; T (:,:,:,end,:,:,:,:,:) += T1(:,:,:,end,:,:,:,:,:) ;
  T1 = T (:,:,1:end-1,:,:,:,:,:,:) ; T1(:,:,end,:,:,:,:,:,:) += T (:,:,end,:,:,:,:,:,:) ;
  T  = T1(:,1:end-1,:,:,:,:,:,:,:) ; T (:,end,:,:,:,:,:,:,:) += T1(:,end,:,:,:,:,:,:,:) ;
  K  = T (1:end-1,:,:,:,:,:,:,:,:) ; K (end,:,:,:,:,:,:,:,:) += T (end,:,:,:,:,:,:,:,:) ;
endfunction

pp = rand(4,4,4);
K = tensor9_opt(pp) ;