MATLAB:如何向量乘两个矩阵数组？

Question

我有两个三维数组,前两个维度表示矩阵,最后一个维度通过参数空间计算,作为一个简单的例子拿

A = repmat([1,2; 3,4], [1 1 4]);

(但假设A(:,:,j)每个人都不同j).一个如何能够容易地进行以每个j两个这样的矩阵阵列的矩阵乘法A和B？

C = A; % pre-allocate, nan(size(A,1), size(B,2)) would be better but slower
for jj = 1:size(A, 3)
  C(:,:,jj) = A(:,:,jj) * B(:,:,jj);
end

当然可以胜任,但是如果第三维更像是1e3元素,那么这非常慢,因为它不使用MATLAB的矢量化.那么,有更快的方法吗？

Answer 1

我现在做了一些时序测试,2x2xN的最快方法是计算矩阵元素:

C = A;
C(1,1,:) = A(1,1,:).*B(1,1,:) + A(1,2,:).*B(2,1,:);
C(1,2,:) = A(1,1,:).*B(1,2,:) + A(1,2,:).*B(2,2,:);
C(2,1,:) = A(2,1,:).*B(1,1,:) + A(2,2,:).*B(2,1,:);
C(2,2,:) = A(2,1,:).*B(1,2,:) + A(2,2,:).*B(2,2,:);

在一般情况下,事实证明for循环实际上是最快的(不要忘记预先分配C!).

如果已经将结果作为矩阵的单元格数组,那么使用cellfun是最快的选择,它也比循环遍历单元格元素更快:

C = cellfun(@mtimes, A, B, 'UniformOutput', false);

但是,必须先调用num2cell first(Ac = num2cell(A, [1 2]))和cell2mat3d数组的情况会浪费太多时间.

这是我为随机的2 x 2 x 1e4设置的时间:

 array-for: 0.057112
 arrayfun : 0.14206
 num2cell : 0.079468
 cell-for : 0.033173
 cellfun  : 0.025223
 cell2mat : 0.010213
 explicit : 0.0021338

显式是指使用2 x 2矩阵元素的直接计算,请参见下文.对于新的随机数组,结果类似,cellfun如果num2cell之前不需要则最快,并且对2x2xN没有限制.对于在第三维上循环的一般3d数组确实是最快的选择.这是时间码:

n = 2;
m = 2;
l = 1e4;

A = rand(n,m,l);
B = rand(m,n,l);

% naive for-loop:
tic
%Cf = nan(n,n,l);
Cf = A;
for jl = 1:l
    Cf(:,:,jl) = A(:,:,jl) * B(:,:,jl);
end;
disp([' array-for: ' num2str(toc)]);

% using arrayfun:
tic
Ca = arrayfun(@(k) A(:,:,k)*B(:,:,k), 1:size(A,3), 'UniformOutput',false);
Ca = cat(3,Ca{:});
disp([' arrayfun : ' num2str(toc)]);

tic
Ac = num2cell(A, [1 2]);
Bc = num2cell(B, [1 2]);
disp([' num2cell : ' num2str(toc)]);

% cell for-loop:
tic
Cfc = Ac;
for jl = 1:l
    Cfc{jl} = Ac{jl} * Bc{jl};
end;
disp([' cell-for : ' num2str(toc)]);

% using cellfun:
tic
Cc = cellfun(@mtimes, Ac, Bc, 'UniformOutput', false);
disp([' cellfun  : ' num2str(toc)]);

tic
Cc = cell2mat(Cc);
disp([' cell2mat : ' num2str(toc)]);

tic
Cm = A;
Cm(1,1,:) = A(1,1,:).*B(1,1,:) + A(1,2,:).*B(2,1,:);
Cm(1,2,:) = A(1,1,:).*B(1,2,:) + A(1,2,:).*B(2,2,:);
Cm(2,1,:) = A(2,1,:).*B(1,1,:) + A(2,2,:).*B(2,1,:);
Cm(2,2,:) = A(2,1,:).*B(1,2,:) + A(2,2,:).*B(2,2,:);
disp([' explicit : ' num2str(toc)]);

disp(' ');

Answer 2

我强烈推荐你使用matlab的MMX工具箱。它可以尽可能快地乘以n维矩阵。

MMX的优点是：

1. 它很容易使用。
2. n维矩阵相乘（实际上可以是二维矩阵数组相乘）
3. 它执行其他矩阵运算（转置、二次乘法、Chol 分解等）
4. 它使用C编译器和多线程计算来加速。

对于这个问题，你只需要编写这样的命令：

C=mmx('mul',A,B);

我在@Amro的答案中添加了以下功能

%# mmx toolbox
function C=func6(A,B,n,m,p)
    C=mmx('mul',A,B);
end

我得到这个结果n=2,m=2,p=1e5：

    1.6571 # FOR-loop
    4.3110 # ARRAYFUN
    3.3731 # NUM2CELL/FOR-loop/CELL2MAT
    2.9820 # NUM2CELL/CELLFUN/CELL2MAT
    0.0244 # Loop Unrolling
    0.0221 # MMX toolbox  <===================

我使用@Amro 的代码来运行基准测试。