MATLAB Answers

3D/2D matrix multiplication without using a loop

2 views (last 30 days)
Pascal Wielsch
Pascal Wielsch on 17 Jun 2021
Commented: Matt J on 18 Jun 2021
Hello dear MATLAB community,
again I have a question to improve the efficiency of my code, by getting rid of the use of loops. Since my programming skills come from using LabVIEW, I often have a harder time using matrix operations instead of loops.
Let me first explain you what I want to do. I have a 3D data matrix "A(i,j,k)", that I want to multiply together with a sensitivity matrix "S" (values stay equal) to get forces and do a conversion from polar to cartesian. I want to do this calculation/conversion for every sub-matrix "A(:,:,k)" to calculate the single forces in the cartesian system. Afterwards I want to calculate the total resulting forces of all sub-matrices and convert them back in a polar system.
The matrices and it's dimension can differ:
  • S_r & S_phi are not always square matrices
  • Row nr. of A is always equal to col nr. of S_r/S_phi
  • Col nr. of A is always equal to length of phipoints
% Input data matrix (3D): A(i,j,k)
A(:,:,1) = [-0.0468 -0.0652 0.2746 0.5272 0.2373 -0.2142 -0.5973 -0.1169
0.0633 0.6673 0.6416 0.3491 -0.0652 -0.0025 -0.8464 -0.7607
-0.0026 0.0000 -0.0373 -0.4378 -0.1612 0.2293 0.2458 0.1656
0.3124 0.2472 0.0048 -0.5034 -0.5051 -0.0434 0.1734 0.3376
-0.0815 -0.4423 0.0308 0.1411 0.1260 0.1881 0.0696 -0.0328];
A(:,:,2) = [-0.0468 -0.0652 0.2746 0.5272 0.2373 -0.2142 -0.5973 -0.1169
0.0633 0.6673 0.6416 0.3491 -0.0652 -0.0025 -0.8464 -0.7607
-0.0026 0.0000 -0.0373 -0.4378 -0.1612 0.2293 0.2458 0.1656
0.0696 0.1881 0.1260 0.1411 0.0308 -0.4423 -0.0815 -0.0328
0.1734 -0.0434 -0.5051 -0.5034 0.0048 0.2472 0.3124 0.3376];
A(:,:,3) = [-0.0468 -0.0652 0.2746 0.5272 0.2373 -0.2142 -0.5973 -0.1169
0.0633 0.6673 0.6416 0.3491 -0.0652 -0.0025 -0.8464 -0.7607
0.1734 -0.0434 -0.5051 -0.5034 0.0048 0.2472 0.3124 0.3376
0.2458 0.2293 -0.1612 -0.4378 -0.0373 0.0000 -0.0026 0.1656
-0.0815 -0.4423 0.0308 0.1411 0.1260 0.1881 0.0696 -0.0328];
% points of all to be calculated segments (index j)
phipoints = 0:2*pi/8:2*pi-2*pi/8;
% Sensitivity matrix: amplitude and phase
S_r = [ 0.0317 0.0378 0.0344 0.0394 0.0333;...
0.0331 0.0410 0.0380 0.0433 0.0375;...
0.0326 0.0415 0.0388 0.0443 0.0393;...
0.0296 0.0386 0.0360 0.0417 0.0383;...
0.0247 0.0334 0.0307 0.0366 0.0354];
S_phi = [ 0 0 0 0 0;...
0 0 0 0 0;...
0 0 0 0 0;...
0 0 0 0 0;...
0 0 0 0 0];
% To be calculate equation for every single sub-matrix A(:,:,k); only for
% documentation; code must be executed using loops
% pol2cart(S_phi+phipoints,S_r*A(:,:,k))
I can calculate this using loops, but my problem is that I have to do this calculation for many sub-matrices (dimension "k" can be up to millions). My current version (see below) takes way to much time, so I'm looking for an alternative way.
%% Current calculation using loop
% Calculation: single forces in cartesian system (X & Y)
for k = 1:size(A,3)
for i = 1:size(A,1)
[Fx(:,:,i,k),Fy(:,:,i,k)] = pol2cart(S_phi(:,i)+phipoints,S_r(:,i)*A(i,:,k));
end
end
% Calculation: resulting forces in X & Y
% - 1st: sum of single forces from all rows (2nd index) of Fx/Fy
% - 2nd: sum of forces of 3rd index of Fx/Fy
Fx_res = sum(sum(Fx,2),3);
Fy_res = sum(sum(Fy,2),3);
% Calculation: total resulting force
[F_res_phi,F_res] = cart2pol(Fx_res,Fy_res); % Conversion: cartesian to polar
Thanks in advance and best regards,
Pascal
  5 Comments
J. Alex Lee
J. Alex Lee on 18 Jun 2021
It's not even clear what "S_phi + phipoints" is supposed to mean whenever S_phi [=] 5xL where L~=5.

Sign in to comment.

Accepted Answer

Pascal Wielsch
Pascal Wielsch on 18 Jun 2021
I finally found a solution to my problem to get rid of the calculation with a loop. You can use the initialization of the input variables from above to test it.
% Extract function variables
I = size(A,1);
J = size(A,2);
K = size(A,3);
% re-arrangement of matrices for calculation (4D)
S_phi_ForCalc = repmat(reshape(S_phi,5,1,I),1,1,1,K);
S_r_ForCalc = repmat(reshape(S_r,5,1,I),1,1,1,K);
A_ForCalc = reshape(permute(A,[2 1 3]),1,J,I,K);
% target equation without using a loop
[Fx,Fy] = pol2cart(S_phi_ForCalc+phipoints,S_r_ForCalc.*A_ForCalc);
% Calculation: resulting forces in X & Y
% - 1st: sum of single forces from all rows (2nd index) of Fx/Fy
% - 2nd: sum of forces of 3rd index of Fx/Fy
Fx_res = sum(sum(Fx,2),3);
Fy_res = sum(sum(Fy,2),3);
% Calculation: total resulting force
[F_res_phi,F_res] = cart2pol(Fx_res,Fy_res); % Conversion: cartesian to polar
  4 Comments
Matt J
Matt J on 18 Jun 2021
Note that in recent Matlab, the repmat's should really not be needed:
S_phi_ForCalc = reshape(S_phi,5,1,I);
S_r_ForCalc = reshape(S_r,5,1,I);

Sign in to comment.

More Answers (2)

J. Alex Lee
J. Alex Lee on 17 Jun 2021
I found this: pagemtimes()
s = rand(5,5);
x = rand(5,8,100000);
tic
y = pagemtimes(s,x);
toc
Elapsed time is 0.019204 seconds.
tic
z = nan(size(x));
for k = 1:size(x,3)
z(:,:,k) = s*x(:,:,k);
end
toc
Elapsed time is 0.261254 seconds.
err = sum(abs(y-z),'all')
err = 1.5855e-10
  2 Comments

Sign in to comment.


Matt J
Matt J on 17 Jun 2021
Edited: Matt J on 17 Jun 2021
Assuming S_r will always be square,
B=reshape( S_r*A(:,:), size(A)); %B(:,:,k)=S_r*A(:,:,k)
  2 Comments
J. Alex Lee
J. Alex Lee on 18 Jun 2021
You would just need to re-figure the size, but this should still work, it's neat!
clc;
clear;
L = 6;
N = 100000;
S_r = rand(5,L);
A = rand(L,8,N);
tic
x = pagemtimes(S_r,A);
toc
Elapsed time is 0.022069 seconds.
tic
y = reshape(S_r*A(:,:),[5,8,N]);
toc
Elapsed time is 0.048928 seconds.
tic
z = nan(5,8,N);
for k = 1:size(A,3)
z(:,:,k) = S_r*A(:,:,k);
end
toc
Elapsed time is 0.237041 seconds.
err = sum(abs(x-y),'all')
err = 0
err = sum(abs(y-z),'all')
err = 1.8814e-10

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!