"Partial" matrix multiplication

9 次查看（过去 30 天）
Hieu Pham
Hieu Pham 2015-6-1
回答： Joss Knight 2015-6-26
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.

请先登录，再回答此问题。

采纳的回答

Azzi Abdelmalek
Azzi Abdelmalek 2015-6-1
编辑：Azzi Abdelmalek 2015-6-1
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))
  1 个评论
Hieu Pham
Hieu Pham 2015-6-2
Thank you for your answer. However, just recently I did a comparison.
function test(D,n)
A = rand(D,3*n);
B = rand(D,3*n);
tic; x = test1(A,B); toc
tic; x = test2(A,B); toc
end
function [out] = test1(A,B)
idx=1:3:size(A,2);
out=cell2mat(arrayfun(@(x) A(:,x:x+2)*B(:,x:x+2)',idx,'un',0));
end
function [out] = test2(A,B)
D = size(B,1);
n = size(B,2) / 3;
out = zeros(D,D*n);
for i = 1:n
out(:,D*(i-1)+1:D*i) = A(:,3*i-2:3*i) * B(:,3*i-2:3*i)';
end
end
And it turns out that your method is (much) slower than a simple for loop...
>> test(1000, 128)
Elapsed time is 5.499530 seconds.
Elapsed time is 1.431573 seconds.

请先登录，再进行评论。

更多回答（2 个）

James Tursa
James Tursa 2015-6-2
编辑：James Tursa 2015-6-2
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:
http://www.mathworks.com/matlabcentral/fileexchange/25977-mtimesx-fast-matrix-multiply-with-multi-dimensional-support
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]));
http://www.mathworks.com/matlabcentral/fileexchange/37515-mmx-multithreaded-matrix-operations-on-n-d-matrices
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:
http://www.mathworks.com/matlabcentral/fileexchange/8773-multiple-matrix-multiplications--with-array-expansion-enabled
  2 个评论
Hieu Pham
Hieu Pham 2015-6-2
编辑：Hieu Pham 2015-6-2
Dear James,
Thank you for your answer. I do have a C compiler so technically I can try mtimesx and mmx. However, at the end of the day, I will have to run my code on a GPU. Are you aware of a way to use mtimesx and/or mmx with GPU? For e.g., a CUDA version?
Last but not least, the last method you proposed seems to be even slower than the one proposed by @Azzi above.
>> test(1000, 128)
Elapsed time is 2.536413 seconds.
Elapsed time is 1.380087 seconds.
Elapsed time is 13.902408 seconds.
The three runtimes correspond to Azzi's method, for loop and your method.
James Tursa
James Tursa 2015-6-2
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
Joss Knight 2015-6-26
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.

类别

在 Help Center 和 File Exchange 中查找有关 Linear Algebra 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by