3D/2D matrix multiplication without using a loop

33 views (last 30 days)
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
Pascal Wielsch
Pascal Wielsch on 18 Jun 2021
I finally found a solution. I will post it in a separate answer.
Thanks a lot for your support.
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
J. Alex Lee
J. Alex Lee on 18 Jun 2021
Well, actually no, the results used to be 5x8xN (not Lx8xN), regardless of L. It doesn't really make sense to me what you mean by S_phi+phipoints when L~=5.
But anyway glad if you've got what you want in the end.
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
Pascal Wielsch
Pascal Wielsch on 17 Jun 2021
Unfortunately I can't find this function. Is it from a toolbox or maybe a later MATLAB version? I use R2019b.
Additionally, I think this will not solve the whole problem, since I still need to calculate this part of the equation: "S_phi(:,i)+phipoints".

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
Pascal Wielsch
Pascal Wielsch on 18 Jun 2021
Thanks. It works for this case, but unfortunately S_r isn't always square. I will edit my question and give a note on the possible dimension of the matrices.
Also, as I mentioned in my recent comment to J. Alex Lee's hint, I still have no clue how to handle the equation "S_phi(:,i)+phipoints" and get rid of the loop.
J. Alex Lee
J. Alex Lee on 18 Jun 2021
Edited: 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.

Categories

Find more on Loops and Conditional Statements in Help Center and File Exchange

Products


Release

R2019b

Community Treasure Hunt

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

Start Hunting!