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,

enter image description here

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

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.

  • $\begingroup$ 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. $\endgroup$ Jul 21 '17 at 17:43
  • $\begingroup$ @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? $\endgroup$ Jul 21 '17 at 17:47
  • 1
    $\begingroup$ 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. $\endgroup$ Jul 21 '17 at 18:08
  • $\begingroup$ 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. $\endgroup$
    – knl
    Jul 23 '17 at 11:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.