如何避免用于矩阵所有排列相乘的for循环？

Question

我有以下代码：

  N=8;
  K=10;
  a=zeros(1,N^(K-1));
  b=zeros(1,N^(K-1));

  for ii=1:K
    p0{ii}=rand(1,N);
    p1{ii}=rand(1,N);
  end

  k=1;
  for j1=1:N
    for j3=1:N
      for j4=1:N
        for j5=1:N
          for j6=1:N
            for j7=1:N
              for j8=1:N
                for j9=1:N
                  for j10=1:N
                    a(k)=p0{1}(j1)*p0{3}(j3)*p0{4}(j4)*p0{5}(j5)*p0{6}(j6)*p0{7}(j7)*p0{8}(j8)*p0{9}(j9)*p0{10}(j10);
                    b(k)=p1{1}(j1)*p1{3}(j3)*p1{4}(j4)*p1{5}(j5)*p1{6}(j6)*p1{7}(j7)*p1{8}(j8)*p1{9}(j9)*p1{10}(j10);
                    k=k+1;
                  end
                end
              end
            end
          end
        end
      end
    end
    end

我无法为N=8 评估此代码，因为它需要很多时间。 p0 和 p1 是大小为 KxN 的矩阵。嵌套的for循环省略了p0和p1的一行，这里第二行对应索引j2。其余的矩阵元素相互相乘。所以总共有N^(K-1) 乘法以获得向量a 和b。

有没有办法在不使用 for 循环的情况下或至少在某个合理的时间内完成这件事？

Answer 1

基本上，您只需将每个p0（或p1）单元格中的每个元素相乘即可。使用reshape 和element-wise multiplication 的一些魔法，这可以简化为一个循环。

让我们看看下面的代码：

N = 3;
K = 10;

for ii = 1:K
  p0{ii} = rand(1, N);
  p1{ii} = rand(1, N);
end  

a = zeros(1, N^(K-1));
b = zeros(1, N^(K-1));

%for ii = 1:K
%  p0{ii} = rand(1, randi(N));
%  p1{ii} = rand(1, randi(N));
%end

tic;
k = 1;
for j1 = 1:N
  for j3 = 1:N
    for j4 = 1:N
      for j5 = 1:N
        for j6 = 1:N
          for j7 = 1:N
            for j8 = 1:N
              for j9 = 1:N
                for j10 = 1:N
                  a(k) = p0{1}(j1)*p0{3}(j3)*p0{4}(j4)*p0{5}(j5)*p0{6}(j6)*p0{7}(j7)*p0{8}(j8)*p0{9}(j9)*p0{10}(j10);
                  b(k) = p1{1}(j1)*p1{3}(j3)*p1{4}(j4)*p1{5}(j5)*p1{6}(j6)*p1{7}(j7)*p1{8}(j8)*p1{9}(j9)*p1{10}(j10);
                  k = k+1;
                end
              end
            end
          end
        end
      end
    end
  end
end
toc;

tic;
aa = p0{1};
bb = p1{1};
% For MATLAB versions R2016 and newer:
for jj = 3:K
  aa = reshape(aa .* p0{jj}.', 1, numel(aa) .* numel(p0{jj}));
  bb = reshape(bb .* p1{jj}.', 1, numel(bb) .* numel(p1{jj}));
end
% For MATLAB versions before R2016b: 
%for jj = 3:K
%  aa = reshape(bsxfun(@times, aa, p0{jj}.'), 1, numel(aa) .* numel(p0{jj}));
%  bb = reshape(bsxfun(@times, bb, p1{jj}.'), 1, numel(bb) .* numel(p1{jj}));
%end
toc;

numel(find(aa ~= a))
numel(find(bb ~= b))

输出：

Elapsed time is 2.39744 seconds.
Elapsed time is 0.00070405 seconds.
ans = 0
ans = 0

看来，a 和 aa 以及 b 和 bb 实际上是相等的，并且建议的解决方案要快得多。我测试了N = 8 只是为了我的解决方案：

Elapsed time is 1.54249 seconds.

如果您通过取消注释相应的行来替换p0 和p1 的初始化，您会看到我的解决方案还允许每个p0（或p1）单元格的长度不同。注意：由于硬编码，这不适用于您的初始解决方案，因此无法在此处进行比较。

另外，请注意，jj = 3:N 在这里也是硬编码的。如果需要省略其他部分，则必须相应修改！

希望有帮助！