3D matrix, multiplication on last dimensions

Question

I have a 3D-matrix A, with size lets say 3x12x100. The first two dimensions define 3×12 matrices, the latter one is simply the linear index. I want a very simple operation on these 100 matrices. For all these matrices, i want them multiplied with its conjugate transpose. With a very simple for loop, i can create this:

data = data;
A = zeros(100, 12, 12);
for i=1:100
    A(i, :, :) = data(:, :, i)'*data(:, :, i);
end

But i like clean code, so i dont really prefer this for-loop. I have done some searching and sometimes find something like mtimesx (which is a custom made MATLAB function from 2010). I think i am missing something very obvious (as usual), because this seems a fairly easy operation (its just an "element-wise" matrix multiplication).

I think your code looks clear and simple and that any other thing that hides the loop will just make it less ovbious — Ander Biguri
Okay, thank you :) Was just wondering, it is such a simple operation (apart from the conjugate transpose). — user2501247
I think you can do this without the loop by appropriately permuteing the matrices and adding singleton dimensions, then applying the dot operator along the right dimension. I don’t think it is worth while to figure out how exactly to write this. The loop code you have will be clearer, and likely faster. — Cris Luengo
My code is however terribly slow. The total number of entries can easily be 800000 (so the result is that many 24x24 matrices), which can easily take around 10 minutes. I have a lot more (simple multiply/add) operations preceding this calculation and these take less than a second. Got any tips? — user2501247

user2501247 user2501247 · Accepted Answer · 2020-04-08T15:22:56

The size of my actual matrix is 3x12x862400. My original script takes about 10 minutes or longer, a variant on what @FangQ posts fixes it in a matter of seconds. My new code is as following, note that it still is under construction and i still need to validate it:

data = rand(3, 12, 862400) + i*rand(3, 12, 862400)
data2 = conj(permute(data, [2 1 3])); % conjugate transpose each matrix
% my data matrix contains 862400 3x12 matrices with complex numbers

Ap = permute(data2, [2 1 4 3]);
Bp = permute(data, [1 4 2 3]);
M = Ap.*Bp;
M = sum(M, 1);
M = permute(M, [2 3 4 1]);

3D matrix, multiplication on last dimensions

2 Answers