26

假设我有一个AxBxC矩阵X 和一个BxD矩阵Y。

是否有一种非循环方法可以将每个C AxB矩阵与相乘Y？

4

10 回答 10

17

作为个人喜好，我希望我的代码尽可能简洁易读。

这是我会做的，尽管它不符合您的“无循环”要求：

for m = 1:C

    Z(:,:,m) = X(:,:,m)*Y;

end

这导致A x D x C矩阵Z。

当然，您始终可以预先分配 Z 以使用Z = zeros(A,D,C);.

于 2009-11-17T05:28:46.653 回答
15

您可以在一行中使用函数NUM2CELL将矩阵X分解为一个单元格数组，并使用 CELLFUN跨单元格进行操作：

Z = cellfun(@(x) x*Y,num2cell(X,[1 2]),'UniformOutput',false);

结果Z是一个1×C元胞数组，其中每个元胞都包含一个A×D矩阵。如果要Z成为A-by-D-by-C矩阵，可以使用CAT函数：

Z = cat(3,Z{:});


注意：我的旧解决方案使用MAT2CELL而不是NUM2CELL，这并不简洁：

[A,B,C] = size(X);
Z = cellfun(@(x) x*Y,mat2cell(X,A,B,ones(1,C)),'UniformOutput',false);
于 2009-11-17T03:34:59.457 回答
8

这是一个单行解决方案（如果要拆分为第三维，则为两个）：

A = 2;
B = 3;
C = 4;
D = 5;

X = rand(A,B,C);
Y = rand(B,D);

%# calculate result in one big matrix
Z = reshape(reshape(permute(X, [2 1 3]), [A B*C]), [B A*C])' * Y;

%'# split into third dimension
Z = permute(reshape(Z',[D A C]),[2 1 3]);

因此现在：Z(:,:,i)包含结果X(:,:,i) * Y

解释：

以上可能看起来令人困惑，但想法很简单。首先，我首先获取的第三个维度X并沿第一个暗淡进行垂直连接：

XX = cat(1, X(:,:,1), X(:,:,2), ..., X(:,:,C))

...困难在于它C是一个变量，因此您无法使用cat或vertcat概括该表达式。接下来我们乘以Y：

ZZ = XX * Y;

最后，我将其拆分回第三维：

Z(:,:,1) = ZZ(1:2, :);
Z(:,:,2) = ZZ(3:4, :);
Z(:,:,3) = ZZ(5:6, :);
Z(:,:,4) = ZZ(7:8, :);

所以你可以看到它只需要一次矩阵乘法，但你必须在前后重塑矩阵。

于 2009-11-16T23:35:10.897 回答
5

我正在处理完全相同的问题，着眼于最有效的方法。我看到了大约三种方法，没有使用外部库（即mtimesx）：

  1. 循环遍历 3D 矩阵的切片
  2. repmat-and-permute 魔法
  3. 细胞乐趣乘法

我最近比较了所有三种方法，看看哪种方法最快。我的直觉是（2）会是赢家。这是代码：

% generate data
A = 20;
B = 30;
C = 40;
D = 50;

X = rand(A,B,C);
Y = rand(B,D);

% ------ Approach 1: Loop (via @Zaid)
tic
Z1 = zeros(A,D,C);
for m = 1:C
    Z1(:,:,m) = X(:,:,m)*Y;
end
toc

% ------ Approach 2: Reshape+Permute (via @Amro)
tic
Z2 = reshape(reshape(permute(X, [2 1 3]), [A B*C]), [B A*C])' * Y;
Z2 = permute(reshape(Z2',[D A C]),[2 1 3]);
toc


% ------ Approach 3: cellfun (via @gnovice)
tic
Z3 = cellfun(@(x) x*Y,num2cell(X,[1 2]),'UniformOutput',false);
Z3 = cat(3,Z3{:});
toc

所有三种方法都产生了相同的输出（呸！），但令人惊讶的是，循环是最快的：

Elapsed time is 0.000418 seconds.
Elapsed time is 0.000887 seconds.
Elapsed time is 0.001841 seconds.

请注意，从一次试验到另一次试验的时间可能会有很大差异，有时 (2) 的结果最慢。这些差异随着更大的数据而变得更加显着。但是对于更大的数据，（3）胜过（2）。循环方法仍然是最好的。

% pretty big data...
A = 200;
B = 300;
C = 400;
D = 500;
Elapsed time is 0.373831 seconds.
Elapsed time is 0.638041 seconds.
Elapsed time is 0.724581 seconds.

% even bigger....
A = 200;
B = 200;
C = 400;
D = 5000;
Elapsed time is 4.314076 seconds.
Elapsed time is 11.553289 seconds.
Elapsed time is 5.233725 seconds.

但是如果循环尺寸比其他方法大得多，则循环方法可能比（2）慢。

A = 2;
B = 3;
C = 400000;
D = 5;
Elapsed time is 0.780933 seconds.
Elapsed time is 0.073189 seconds.
Elapsed time is 2.590697 seconds.

因此，在这种（可能是极端的）情况下，（2）以一个很大的因素获胜。可能没有一种方法在所有情况下都是最佳的，但循环仍然非常好，并且在许多情况下是最好的。它在可读性方面也是最好的。绕开！

于 2015-10-21T21:18:14.790 回答
1

没有。有几种方法，但它总是直接或间接地循环出现。

只是为了满足我的好奇心，你为什么要这样呢？

于 2009-11-16T22:49:36.597 回答
1

要回答这个问题并为了便于阅读，请参阅：

  • ndmult，作者 ajuanpi (Juan Pablo Carbajal)，2013，GNU GPL

输入

  • 2 个阵列
  • 暗淡

例子

 nT = 100;
 t = 2*pi*linspace (0,1,nT)’;

 # 2 experiments measuring 3 signals at nT timestamps
 signals = zeros(nT,3,2);
 signals(:,:,1) = [sin(2*t) cos(2*t) sin(4*t).^2];
 signals(:,:,2) = [sin(2*t+pi/4) cos(2*t+pi/4) sin(4*t+pi/6).^2];

 sT(:,:,1) = signals(:,:,1)’;
 sT(:,:,2) = signals(:,:,2)’;
   G = ndmult (signals,sT,[1 2]);

来源

原始来源。我添加了内联注释。

function M = ndmult (A,B,dim)
  dA = dim(1);
  dB = dim(2);

  # reshape A into 2d
  sA = size (A);
  nA = length (sA);
  perA = [1:(dA-1) (dA+1):(nA-1) nA dA](1:nA);
  Ap = permute (A, perA);
  Ap = reshape (Ap, prod (sA(perA(1:end-1))), sA(perA(end)));

  # reshape B into 2d
  sB = size (B);
  nB = length (sB);
  perB = [dB 1:(dB-1) (dB+1):(nB-1) nB](1:nB);
  Bp = permute (B, perB);
  Bp = reshape (Bp, sB(perB(1)), prod (sB(perB(2:end))));

  # multiply
  M = Ap * Bp;

  # reshape back to original format
  s = [sA(perA(1:end-1)) sB(perB(2:end))];
  M = squeeze (reshape (M, s));
endfunction
于 2014-01-08T23:45:30.513 回答
1

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

MMX的优点是：

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

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

C=mmx('mul',X,Y);

这是所有可能方法的基准。有关更多详细信息，请参阅此问题。

    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  <===================
于 2017-08-12T06:06:03.627 回答
0

我认为递归，但这是您可以做的唯一其他非循环方法

于 2009-11-16T22:39:13.420 回答
0

您可以“展开”循环，即按顺序写出循环中发生的所有乘法

于 2009-11-16T22:42:24.137 回答
0

我想分享我对以下问题的回答：

1）制作两个张量（任何价）的张量积；

2) 使两个张量沿任意维度收缩。

这是我的第一个和第二个任务的子程序：

1）张量积：

function [C] = tensor(A,B)
   C = squeeze( reshape( repmat(A(:), 1, numel(B)).*B(:).' , [size(A),size(B)] ) );
end

2) 收缩：这里 A 和 B 是分别沿着维度 i 和 j 收缩的张量。当然，这些尺寸的长度应该相等。没有对此进行检查（这会掩盖代码），但除此之外它运行良好。

   function [C] = tensorcontraction(A,B, i,j)
      sa = size(A);
      La = length(sa);
      ia = 1:La;
      ia(i) = [];
      ia = [ia i];

      sb = size(B);
      Lb = length(sb);
      ib = 1:Lb;
      ib(j) = [];
      ib = [j ib];

      % making the i-th dimension the last in A
      A1 = permute(A, ia);
      % making the j-th dimension the first in B
      B1 = permute(B, ib);

      % making both A and B 2D-matrices to make use of the
      % matrix multiplication along the second dimension of A
      % and the first dimension of B
      A2 = reshape(A1, [],sa(i));
      B2 = reshape(B1, sb(j),[]);

      % here's the implicit implication that sa(i) == sb(j),
      % otherwise - crash
      C2 = A2*B2;

      % back to the original shape with the exception
      % of dimensions along which we've just contracted
      sa(i) = [];
      sb(j) = [];
      C = squeeze( reshape( C2, [sa,sb] ) );
   end

有批评者吗？

于 2019-08-25T17:04:28.047 回答