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

1 view (last 30 days)
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?

Accepted Answer

DGM
DGM on 23 Oct 2021
Edited: DGM on 23 Oct 2021
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.
  2 Comments
Qiu Xu
Qiu Xu on 23 Oct 2021
Thanks very much for your help. The second method is very fast for general cases.
Qiu Xu
Qiu Xu on 23 Oct 2021
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.

Sign in to comment.

More Answers (1)

Chunru
Chunru on 23 Oct 2021
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

Categories

Find more on Performance and Memory in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!