each columun of a 2D matrix multiply with a 3D array

Hi,

I have a 2D matrix and a 3D array, I want to perform the following calculations:

 A=rand(200,300);
   B=rand(60,60,200);
   c=0;
   for j=1:300
        d=0;
        for k=1:200
              d=d+A(k,j)*B(:,:,k);
        end
        c=c+d;
   end

The a two for loops, the calculations are very slowly. Are there some good methods to aviod the loops?

Answer 1

[DGM]

This is something, but I bet there are other ways too. This is about 10x faster on my hardware/version, but only about 2-3x faster on the web version.

% array size parameters
pagesz = [10 10];
npages = 10;
wa = 10;

% build test arrays
A = rand(npages,wa);
B = rand(pagesz(1),pagesz(2),npages);

% METHOD 1
tic
c = 0;
for j = 1:wa
    d = 0;
    for k = 1:npages
        d = d+A(k,j)*B(:,:,k);
    end
    c = c+d;
end
toc

Elapsed time is 0.006433 seconds.

% METHOD 2
tic
c2 = sum(B.*permute(sum(A,2),[3 2 1]),3);
toc

Elapsed time is 0.004335 seconds.

% error metric (should be approximately zero)
immse(c,c2)

ans = 1.1517e-29

% METHOD 3 (Chunru's proposal)
tic
c3 = 0;
for j = 1:wa
    c3 = c3 + sum(reshape(A(:, j), 1, 1, []).*B, 3);
end
toc

Elapsed time is 0.006630 seconds.

% error metric (should be approximately zero)
immse(c,c3)

ans = 0

While Chunru's proposal might not be generally as fast as what I proposed, it preserves the order of operations, and so the accumulated error will be smaller.

[Qiu Xu]

Thanks very much for your help. The second method is very fast for general cases.

[Qiu Xu]

Dear DGM,

I have another promblem. Now I have two 2D matrix u and v, and a 3D array B. I want to perform the following calculations, for example

u=rand(200,300);

v=rand(200,300);

B=rand(60,60,300)

c=0;

for j=1:300

u1=0;

v1=0;

for k=1:200

u1=u1+u(k,j)*B(:,:,k);

v1=v1+v(k,j)*B(:,:,k);

end

c=c+u1.*v1;

end

I do not know how to aviod two loops due to the term u1.*v1 in c.

Thanks a lot.

Answer 2

[Chunru]

A=rand(20,30);
B=rand(6,6,20);
c=0;
for j=1:30
    d=0;
    for k=1:20
          d=d+A(k,j)*B(:,:,k);
    end
    c=c+d;
end
whos

  Name       Size              Bytes  Class     Attributes

  A         20x30               4800  double              
  B          6x6x20             5760  double              
  c          6x6                 288  double              
  d          6x6                 288  double              
  j          1x1                   8  double              
  k          1x1                   8  double              

c

c = 6×6

142.8301 142.6174 164.7719 150.3854 156.3244 157.9079 169.1241 129.4667 141.7922 163.7245 165.0511 136.1945 141.9052 173.1903 127.6732 109.3965 144.3338 175.1451 125.7923 162.6825 183.5525 150.5645 139.2519 129.0993 134.5237 171.8890 156.9149 146.8089 163.1632 133.5324 162.8775 145.6163 152.5322 160.0175 163.6295 147.2230

c=0;
for j=1:30
    c = c + sum(reshape(A(:, j), 1, 1, []).*B, 3);
end
c

c = 6×6

142.8301 142.6174 164.7719 150.3854 156.3244 157.9079 169.1241 129.4667 141.7922 163.7245 165.0511 136.1945 141.9052 173.1903 127.6732 109.3965 144.3338 175.1451 125.7923 162.6825 183.5525 150.5645 139.2519 129.0993 134.5237 171.8890 156.9149 146.8089 163.1632 133.5324 162.8775 145.6163 152.5322 160.0175 163.6295 147.2230

[Qiu Xu]

Thanks very much