How do I vectorize this 4-D matrix operation in Matlab?

Question

I need to do a large number of matrix operations in Matlab, but one of the matrices is 4-D and my normal instincts for correct vectorization are failing me. Right now I'm using a loop. Maybe somebody can fix my dumb brain?

for k=1:kNum
   for l=1:lNum
      L(k,l,:) = squeeze(M(k,l,:,:))\squeeze(P(k,l,:));
   end
end

Answer 1

I don't think that it is possible to vectorize this expression very much. Generally, vectorization is used with assignments to vectors along with functions which output seperate results for each item, for instance

for i = 1:N
   A(i) = sin(i)
end

is equal to

i = 1:N
A(i) = sin(i)

because sin does elementwise sine of the elements. When solving linear equation sets using mldivide/"backslash operator" the solution depends on all the elements in the vector - you can't just do element wise mldivide unless the matrix has special properties.

For optimizing this loop, you could consider using parfor-loops (the looping variables are normal numbers incremented by one, so it should be as simple as doing "parfor" instead of "for). This requires the parallel computing toolbox. You could also unroll the loops to get better performance, if the numbering of the indices is somewhat predictible.

Or, you could try to write an .mex-file doing the same operation. In C/Fortran looping is as fast as / faster than MATLAB vectorization.

Another thing to consider is reordering the memory access: Since matrices are just long arrays accessed using a [i + j*dim1 + k*dim1*dim2] pattern it is important to ensure that the inner loop corresponds to the memory which is stored sequentially. For 2D matrices, at least, doing

for k=1:kNum
   for l=1:lNum
       L(k,l,:) = squeeze(M(k,l,:,:))\squeeze(P(k,l,:));

is slower than

for l= 1:lNum
   for k=1:kNum
       L(k,l,:) = squeeze(M(k,l,:,:))\squeeze(P(k,l,:));

since the elements corresponding to the k-index are stored sequentially. Mathworks has some information on it, but it is not really something specific to MATLAB per se, any somewhat advanced text on arrays for C or FORTRAN should go into more details. There is also some cache stuff going on under the hood here.

edit: While not useful here as the right hand side and left hand side is different in each solution, if you are doing something like

for i = 1:N
    solution(:,i) = A\b(:,i)
end

you might as well do

solution = A\b