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

Question

I have the following code:

  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

I am not able to evaluate this code for N=8, because it takes a lot of time. p0 and p1 are matrices of size KxN. The nested for loop is omiting one row of p0 and p1, here the second row corresponding to the index j2. The rest of the matrix elements are multiplied with each other. So in total there are N^(K-1) multiplications in order to obtain the vectors a and b.

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

Answer 1

Basically, you just multiply each element from each p0 (or p1) cell with each other. Using some magic from reshape and element-wise multiplication, this can simplified to one loop.

Let's have a look at the following code:

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

The output:

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

It seems, a and aa as well as b and bb are actually equal, and the proposed solution is much faster. I tested N = 8 just for my solution:

Elapsed time is 1.54249 seconds.

If you replace the initialization of p0 and p1 by uncommenting the corresponding lines, you'll see my solution also allows for varying lengths for each p0 (or p1) cell. Notice: This doesn't work for your initial solution due to the hardcoding, so comparisons can't be made here.

Also, please note, that jj = 3:N is hardcoded here, too. If other parts should be omitted instead, this has to be modified accordingly!

Hope that helps!