logo标签列表

MATLAB：如何对两个矩阵数组进行向量乘法？

arraysmatlabmatrixvectorizationmatrix-multiplication
8

我有两个三维数组，前两维代表矩阵，最后一维遍历参数空间。以一个简单的例子说明：

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

假设对于每个j，A(:,:,j)都不同。如何轻松地对这样的矩阵数组A和B执行按j矩阵乘法？

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的矢量化。那么，有更快的方法吗？

你是否实际计时了循环？对于最近的Matlab版本，它可能非常快。你期望“矢量化”版本会快多少？谢谢。 - eat
@eat：对于1000个参数，它是7的因子（MATLAB R2010a），我正在优化循环中使用它，所以这很重要-我现在找到了一个解决方案，午饭后我会发布它。 - Tobias Kienzler
1
可能是将3D矩阵与2D矩阵相乘的重复问题。 - gnovice
@TobiasKienzler：我猜你是在预分配矩阵C吗？ - Amro
5个回答
6

我现在进行了一些计时测试，发现对于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 (Ac = num2cell(A, [1 2]))，再使用 cell2mat，这会浪费太多时间。

这是我为一个随机的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

Explicit是指直接计算2 x 2矩阵元素，见下文。对于新的随机数组，结果是相似的，如果不需要在之前使用num2cell并且没有限制到2x2xN，cellfun是最快的。对于一般的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(' ');
1
确实聪明。你可能以后需要接受自己的答案；）。谢谢 - eat
1
不要被CELLFUN所迷惑，它里面有一个隐藏的循环...因此，更简单的方法是直接编写：C = arrayfun(@(k) A(:,:,k)*B(:,:,k), 1:size(A,3), 'UniformOutput',false); C = cat(3,C{:});。这两种方法都不比原始的for循环更好！ - Amro
@Amro：你是对的，我现在进行了时间测试。 arrayfun 的速度几乎与 num2cell + cellfun + cell2mat 相同，结果原始的 for 循环确实是最快的（是的，我预先分配了 C），除非您已经有单元。 - Tobias Kienzler
1
@TobiasKienzler：我发布了一些自己的基准测试...正如预期的那样，FOR循环非常快，特别是在最近版本的MATLAB中使用即时（JIT）加速器进行改进后。 - Amro
4
这是我的基准测试，比较了@TobiasKienzler答案中提到的方法。我使用TIMEIT函数来获取更准确的时间。
function [t,v] = matrixMultTest()
    n = 2; m = 2; p = 1e5;
    A = rand(n,m,p);
    B = rand(m,n,p);

    %# time functions
    t = zeros(5,1);
    t(1) = timeit( @() func1(A,B,n,m,p) );
    t(2) = timeit( @() func2(A,B,n,m,p) );
    t(3) = timeit( @() func3(A,B,n,m,p) );
    t(4) = timeit( @() func4(A,B,n,m,p) );
    t(5) = timeit( @() func5(A,B,n,m,p) );

    %# check the results
    v = cell(5,1);
    v{1} = func1(A,B,n,m,p);
    v{2} = func2(A,B,n,m,p);
    v{3} = func3(A,B,n,m,p);
    v{4} = func4(A,B,n,m,p);
    v{5} = func5(A,B,n,m,p);
    assert( isequal(v{:}) )
end

%# simple FOR-loop
function C = func1(A,B,n,m,p)
    C = zeros(n,n,p);
    for k=1:p
        C(:,:,k) = A(:,:,k) * B(:,:,k);
    end
end

%# ARRAYFUN
function C = func2(A,B,n,m,p)
    C = arrayfun(@(k) A(:,:,k)*B(:,:,k), 1:p, 'UniformOutput',false);
    C = cat(3, C{:});
end

%# NUM2CELL/FOR-loop/CELL2MAT
function C = func3(A,B,n,m,p)
    Ac = num2cell(A, [1 2]);
    Bc = num2cell(B, [1 2]);
    C = cell(1,1,p);
    for k=1:p
        C{k} = Ac{k} * Bc{k};
    end;
    C = cell2mat(C);
end

%# NUM2CELL/CELLFUN/CELL2MAT
function C = func4(A,B,n,m,p)
    Ac = num2cell(A, [1 2]);
    Bc = num2cell(B, [1 2]);
    C = cellfun(@mtimes, Ac, Bc, 'UniformOutput', false);
    C = cell2mat(C);
end

%# Loop Unrolling
function C = func5(A,B,n,m,p)
    C = zeros(n,n,p);
    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,:);
end

结果为：

结果：

>> [t,v] = matrixMultTest();
>> t
t =
      0.63633      # FOR-loop
      1.5902       # ARRAYFUN
      1.1257       # NUM2CELL/FOR-loop/CELL2MAT
      1.0759       # NUM2CELL/CELLFUN/CELL2MAT
      0.05712      # Loop Unrolling

如我在评论中所解释的那样，一个简单的FOR循环是最好的解决方案（除了在最后一种情况下进行循环展开，这仅适用于这些小的2x2矩阵）。

恐怕你刚被Ali的回答抢了勾选标记，他介绍了在2012年之前不存在的MMX工具箱... - Tobias Kienzler
@TobiasKienzler 哦，那没关系。毕竟，很难超越C代码！我查看了MMX工具箱的源代码，它基本上是创建线程（与处理器数量相同），每个线程调用一个矩阵乘法函数来处理其分配的矩阵切片。如果您在编译时启用了优化，则会使用dgemm BLAS例程（从Intel MKL库中提供的MATLAB）来执行矩阵乘法，这是MATLAB内部使用的相同例程。 - Amro
1
对于小型的2x2矩阵而言，你需要注意过度订阅（随MATLAB一起提供的MKL本身是多线程的，同时MMX工具箱也从多个线程调用它）。通过使用针对小矩阵乘法优化的库，你可能会获得更好的性能（BLAS在大矩阵上效果突出）。你可以从Ali的计时中看到这一点；MMX几乎与循环展开版本的运行时间相同。现在想象一下用C实现相同的代码！在我看来，问题是受限于内存而不是CPU，并且线程在这里的效果较差，一切都取决于良好的缓存重用。 - Amro
3
我强烈推荐您使用Matlab的MMX工具箱。它能够以最快的速度相乘n维矩阵。

MMX的优势包括：

  1. 它 易于使用。
  2. 可以相乘n维矩阵（实际上它可以相乘2-D矩阵数组）
  3. 它执行其他矩阵运算（转置、Quadratic Multiply、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的代码来运行基准测试。

1
我的一个旧问题得到了不错的更新 :-) 当然，那个工具箱在2012年之前是不存在的... - Tobias Kienzler
1

根据我的经验，更快的方法是使用三维矩阵的点乘和求和。下面的函数function z_matmultiply(A,B)用于乘以两个深度相同的三维矩阵。点乘尽可能并行地执行，因此您可能需要在大量重复操作时检查此函数的速度并将其与其他函数进行比较。

function C = z_matmultiply(A,B)

[ma,na,oa] = size(A);
[mb,nb,ob] = size(B);

%preallocate the output as we will do a loop soon
C = zeros(ma,nb,oa);

%error message if the dimensions are not appropriate
if na ~= mb || oa ~= ob
    fprintf('\n z_matmultiply warning: Matrix Dimmensions Inconsistent \n')
else

% if statement minimizes for loops by looping the smallest matrix dimension 
if ma > nb
    for j = 1:nb
        Bp(j,:,:) = B(:,j,:);
        C(:,j,:) = sum(A.*repmat(Bp(j,:,:),[ma,1]),2);
    end
else
    for i = 1:ma
        Ap(:,i,:) = A(i,:,:);
        C(i,:,:) = sum(repmat(Ap(:,i,:),[1,nb]).*B,1);
    end 
end

end
1
你可以使用bsxfun代替repmat。 - Shai
2
在Matlab中最好不要使用i和j作为变量名。 - Shai
1
一种技术是创建一个2Nx2N的稀疏矩阵，并在对角线上嵌入2x2矩阵，分别为A和B。使用稀疏矩阵进行乘积运算，并通过巧妙的索引获取结果，然后将其重新整形为2x2xN。
但我怀疑这不会比简单循环更快。
好主意，虽然你的疑虑可能是正确的。如果你感兴趣的话，我找到了一个使用cellfun的解决方案。 - Tobias Kienzler
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接