# Vectorizing Matrix Multiplication

I would like to do the following operation: I have a "4D" matrix A and a "3D" matrix B. Both A and B are actually 2D matrices, where for A, each element is a 2D matrix, and for B, each element is a 1D matrix or vector. I would like to, for each element in the 2D floor plan, multiply the matrix on A and the matrix in B. Graphically,

Currently, my code is:

% size(A) = Nx, Ny, n, n
% size(B) = Nx, Ny, n
% size(C) = Nx, Ny, n

for k1 = 1 : Nx
for k2 = 1 : Ny
C(k1,k2,:) = squeeze(A(k1,k2,:,:)) * squeeze(B(k1,k2,:));
end
end


I am looking for a way to "vectorize" the code and avoid loops completely. I am hoping such a vectorization will speed up my program.

• Vectorizing this kind of computation would lead to a bunch of large temporary arrays which is not very good if you care about the performance. These will be larger than normal since each will be a matrix instead of just a vector. You may want to write this part in C, Fortran, or Julia as non-allocating with the appropriate loops, and link it through MEX if you really need it in MATLAB. – Chris Rackauckas Jul 21 '17 at 17:43
• @ChrisRackauckas For a typical computation, Nx = Ny = 128, n = 6. Would you recommend using MEX? How much of a speed increase could I expect to see? – Geoffrey Xiao Jul 21 '17 at 17:47
• It's hard to know what kind of speedup to expect because it's dependent on the size. Usually you just have to try it to see. But actually, looking at the code again, your biggest problem is that your indexing is backwards. A(k1,k2,:,:). MATLAB is column-major, so that kind of operation is faster if it's A(:,:,k1,k2) since you always want to index down columns (i.e. first indices not latter indices). I'd fix that first and see if that's enough of a speed boost to be okay with. – Chris Rackauckas Jul 21 '17 at 18:08
• I recommend rewriting your code so that you make for loops over $n$ which is the smallest dimension. You should get quite good speedups that way. There is not any out-of-the-box solution for vectorizing n-dimensional arrays as such operations are not intrisic to BLAS. – knl Jul 23 '17 at 11:12