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

[Matt J]

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]

Yes, you're right. The last line of the first code section (k=1; pol2cart(S_phi+phipoints,S_r*A(:,:,k))) should only be for documentation purposes. I've corrected it.

Now it should be possible to run the first section of the code (setting up variables) and afterwards the second section with my current calculation using loops.

[Matt J]

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]

I finally found a solution. I will post it in a separate answer.

Thanks a lot for your support.

[J. Alex Lee]

It's not even clear what "S_phi + phipoints" is supposed to mean whenever S_phi [=] 5xL where L~=5.

Answer 1

[Pascal Fenstermann]

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

[J. Alex Lee]

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]

Yes, but it's basically the same. It's just rearranged for the calculation.

[J. Alex Lee]

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]

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);

Answer 2

[J. Alex Lee]

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

[Pascal Fenstermann]

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".

[J. Alex Lee]

Oops, pagemtimes is available 2020b+.

Answer 3

[Matt J]

Assuming S_r will always be square,

 B=reshape( S_r*A(:,:), size(A));  %B(:,:,k)=S_r*A(:,:,k)

[Pascal Fenstermann]

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]

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