Matrix multiply slices of 3d Matricies

Given two 3d matricies, A and B with

size(A) = (n, m, k)

and

size(B) = (m, p, k)

perform matrix multiplications on each slice obtained by fixing the last index, yielding a matrix C with

size(C) = (n, p, k).

To clarify, we would have

C(:, :, 1) = A(:, :, 1)*B(:, :, 1), ..., C(:, :, k) = A(:, :, k)*B(:, :, k).

I need to do this with gpuArrays in the most efficient manner possible.

Answer 1

[Jill Reese]

If you have MATLAB R2013b, you can use the new gpuArray pagefun function like so:

C = pagefun(@mtimes, A, B);

[Dan Ryan]

Brilliant! Thanks!

Answer 2

[James Tursa]

If you are not restricted to gpuArrays you can do this:

C = mtimesx(A,B);

The MTIMESX function passes pointers to the slice data to BLAS library functions in the background, so it is pretty fast. You can find MTIMESX here:

http://www.mathworks.com/matlabcentral/fileexchange/25977-mtimesx-fast-matrix-multiply-with-multi-dimensional-support

MTIMESX is not yet multi-threaded across the third dimension (but an update is in the works). A nD matrix multiply multi-threaded on the third dimension called MMX can also be used:

C = MMX('mult', A, B);

MMX can be found here:

http://www.mathworks.com/matlabcentral/fileexchange/37515-mmx-multithreaded-matrix-operations-on-n-d-matrices

[Dan Ryan]

great suggestion, I will keep an eye on this project

Answer 3

[Azzi Abdelmalek]

n=3;
m=4;
k=5;
p=2;
A=rand(n,m,k)
B=rand(m,p,k)
C=zeros(n,p,k)
for ii=1:k
 C(:,:,ii)=A(:,:,ii)*B(:,:,ii)
end

[Sean de Wolski]

Could also use a parfor loop since each iteration is independent of others.

[Dan Ryan]

I was hoping to see some fully parallel version that could be implemented in one shot on the gpu. This would mean no looping constructs.

[Azzi Abdelmalek]

I don't see how you could do it without a for loop.

[Sean de Wolski]

What's wrong with a loop?

[Dan Ryan]

Big slowdown... for example:

A = gpuArray.rand(1000, 100, 100, 'single');
B = gpuArray.rand(1000, 100, 'single');
C = gpuArray.zeros(1000, 100, 100, 'single');

Compare

for idx = 1:100
   C(:, :, idx) = A(:, :, idx).*B;
end

with

C = bsxfun(@times, A, B);

There is about a factor of 50 slowdown with the for loop.

[Azzi Abdelmalek]

bsxfun don't work with

n=3;
m=4;
k=5;
p=2;
A=rand(n,m,k)
B=rand(m,p,k)
C = bsxfun(@times, A, B);

[Jill Reese]

Dan,

Can you elaborate on the sizes of m, n , k, and p that you are interested in? It would be useful to know a ballpark number for the size of problem you want to solve. Do you have many small page sizes, a few large pages, or something else?

Thanks,

Jill