I've recently vectorized most of the code for an optimization problem but this one still eludes me. As the procedure and not the actual values are important here, I have replaced the original matrices with random dummies.

Given a matrix A of shape (m x n) and a three dimensional matrix B of shape (n x (n-1)*r x m) multiply the i-th row of A with the i-th slice/page of B for i=1,...,m . This yields a matrix of shape m x (n-1)*r.
Note that I have arranged the m two dimensional matrices in a three dimensional matrix B, which is not necessarily a given. B can as well be modeled as a block matrix or other applicable shapes.

As this code is evaluated very often as part of the optimization, it is crucial that it is as fast as possible. Memorywise it does not need to be optimal, but it should remain within feasible limits, as m=6000, n=6 and r=2 are rather the lower end of the expected dimension of the problem.

What I have tried:

The code below is the first solution that came to my mind and the one that so far has proven the fastest. I have tried replacing the for-loop by using cellfun, but this has not proven to be faster.
Trying to approach this in a vectorized fashion, I also tried structuring B instead of a 3D matrix as a 2D, horizontal block matrix B* with the 'slices' of B simply appended as columns. Then matrix multiplication can then be used to obtain a (m x (n-1)*r*m) matrix where the solution is found on the diagonal + some offset. However, this is of course horribly inefficient, as most of the result matrix isn't required and such a matrix becomes intractable from a memory perspective pretty quickly.

m = 6000;
n = 6;
r = 2;
A = rand(m,n);
B = rand(n,(n-1)*r,m);
result = zeros(m,(n-1)*r);
for i=1:m
    result(i,:) = A(i,:) * B(:,:,i);
  • \$\begingroup\$ cellfun is implemented as an M-file with a loop, so it’s unsurprising that it’s not any faster. I have no comments on your code, it’s correct and likely optimal. Loops are no longer slow in MATLAB, vectorization is only beneficial if it doesn’t involve large intermediate matrices or complex indexing. \$\endgroup\$ Mar 24, 2019 at 0:36
  • \$\begingroup\$ Could it be you meant: result = zeros(m,(n-1)*r, m); and result(:, :, i) = A(i,:) * B(:,:,i);? \$\endgroup\$
    – Royi
    Apr 26, 2019 at 12:43


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.