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

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.

Is there any way to do this thing without using for loops or at least in some reasonable time?

Answer 1

基本上，您只需将每个 p0（或 p1）单元格中的每个元素相互相乘即可。使用 reshape and 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 在这里也是硬编码的。如果要省略其他部分，则必须相应修改！

希望对您有所帮助！

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

How to avoid for loops for multiplication of all permutations of a matrix?

matlab

for-loop

simplification