"Partial" matrix multiplication
7 views (last 30 days)
Show older comments
Suppose that I have two matrices, A and B, both have size Dx(3N). I want to multiply each block of 3 consecutive columns in A with the transpose of the corresponding block of 3 consecutive columns in B (the result of each of these multiplications would be a DxD matrix). What are the best ways to do this?
For example, let's say
A = [a_1, a_2, a_3, b_1, b_2, b_3, c_1, c_2, c_3]
B = [x_1, x_2, x_3, y_1, y_2, y_3, z_1, z_2, z_3]
where a_i, b_i, c_i, x_i, y_i, z_i all have size Dx1. I want to compute
[a_1, a_2, a_3]*[x_1, x_2, x_3]'
[b_1, b_2, b_3]*[y_1, y_2, y_3]'
[c_1, c_2, c_3]*[z_1, z_2, z_3]'
and of course, I need to store the results.
0 Comments
Accepted Answer
Azzi Abdelmalek
on 1 Jun 2015
Edited: Azzi Abdelmalek
on 1 Jun 2015
A=randi(9,3,9)
B=randi(9,3,9)
idx=1:3:size(A,2)
out=cell2mat(arrayfun(@(x) A(:,x:x+2)*B(:,x:x+2)',idx,'un',0))
More Answers (2)
James Tursa
on 2 Jun 2015
Edited: James Tursa
on 2 Jun 2015
If you have a C compiler installed, you can use an FEX submission called mtimesx which does nD matrix multiply with built-in transpose capability (does a virtual transpose, not an actual transpose):
[m,n] = size(A);
n3 = n/3;
Ar = reshape(A,m,3,n3);
Br = reshape(B,m,3,n3);
C = mtimesx(Ar,Br,'t','speedomp');
You can find mtimesx here:
Another option is mmx, but you will have to do the nD transpose manually via a permute:
[m,n] = size(A);
n3 = n/3;
Ar = reshape(A,m,3,n3);
Br = reshape(B,m,3,n3);
C = mmx(Ar,permute(Br,[2 1 3]));
If you don't have a C compiler installed, you can use a different m-file based routine called multiprod:
[m,n] = size(A);
n3 = n/3;
Ar = reshape(A,m,3,n3);
Br = reshape(B,m,3,n3);
C = multiprod(Ar,permute(Br,[2 1 3]));
You can find multiprod here:
2 Comments
James Tursa
on 2 Jun 2015
None of these methods have CUDA versions to my knowledge. For CUDA, you may need to write the code from scratch. If the row size is not too big, hand coding the individual (m x 3) * (m x 3)' multiplies directly element-by-element might be faster than using loops.
Joss Knight
on 26 Jun 2015
If you're running this on a GPU using Parallel Computing Toolbox, as you say, then you can use pagefun:
rows = size(A,1);
assert(size(B,1) == rows);
A = reshape(gpuArray(A), rows, 3, []);
B = reshape(gpuArray(B), rows, 3, []);
Bt = pagefun(@transpose, B);
C = pagefun(@mtimes, A, Bt);
The result C is a rows x rows x (cols/3) ND array.
There is currently no equivalent to pagefun for the CPU, but the CPU will work fine with a loop.
0 Comments
See Also
Categories
Find more on GPU Computing in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!