3D/2D matrix multiplication without using a loop
显示 更早的评论
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 个评论
Matt J
2021-6-18
Your example code doesn't run, as I've now demonstrated above by launching it within your post. Therefore, it's not completely clear what the code is supposed to do.
Pascal Fenstermann
2021-6-18
Matt J
2021-6-18
Yes, it does run, but the operation done in your documentation line
pol2cart(S_phi+phipoints,S_r*A(:,:,k))
does not match what you do later,
for i = 1:size(A,1)
[Fx(:,:,i,k),Fy(:,:,i,k)] = pol2cart(S_phi(:,i)+phipoints,S_r(:,i)*A(i,:,k));
end
Note that S_r(:,i)*A(i,:,k) always produces a rank-1 matrix as input to pol2cart, whereas S_r*A(:,:,k) is a general matrix.
Pascal Fenstermann
2021-6-18
J. Alex Lee
2021-6-18
It's not even clear what "S_phi + phipoints" is supposed to mean whenever S_phi [=] 5xL where L~=5.
采纳的回答
Pascal Fenstermann
2021-6-18
4 个评论
J. Alex Lee
2021-6-18
Didn't your results used to be of same dimensions as A? Now they are 4D...you sure that is what you needed?
Pascal Fenstermann
2021-6-18
J. Alex Lee
2021-6-18
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
2021-6-18
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);
更多回答（2 个）
I found this: pagemtimes()
s = rand(5,5);
x = rand(5,8,100000);
tic
y = pagemtimes(s,x);
toc
tic
z = nan(size(x));
for k = 1:size(x,3)
z(:,:,k) = s*x(:,:,k);
end
toc
err = sum(abs(y-z),'all')
Assuming S_r will always be square,
B=reshape( S_r*A(:,:), size(A)); %B(:,:,k)=S_r*A(:,:,k)
2 个评论
Pascal Fenstermann
2021-6-18
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
tic
y = reshape(S_r*A(:,:),[5,8,N]);
toc
tic
z = nan(5,8,N);
for k = 1:size(A,3)
z(:,:,k) = S_r*A(:,:,k);
end
toc
err = sum(abs(x-y),'all')
err = sum(abs(y-z),'all')
类别
在 帮助中心 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
