3D行列と2D行列を乗算する

https://stackoverflow.com/questions/1745299

20-09-2019
|

質問

私が持っていると仮定します ア×B×C マトリックス X そして BxD マトリックス Y.

それぞれを乗算できる非ループ方法はありますか? C 斧B 行列 Y?

解決

関数を使用してこれを 1 行で実行できます NUM2CELL マトリックスを壊すために X cell 配列に変換し、セルファンセル全体で動作するには:

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

結果 Z です 1 バイ C 各セルに エーバイディー マトリックス。あなたが望むなら Z になる A×D×C マトリックスを使用できます。猫関数：

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

注記： 私が使用していた古いソリューション MAT2セルの代わりに NUM2CELL, 、それほど簡潔ではありませんでした:

[A,B,C] = size(X);
Z = cellfun(@(x) x*Y,mat2cell(X,A,B,ones(1,C)),'UniformOutput',false);

他のヒント

個人的な好みとして、私はコードができるだけ簡潔で読みやすいものであることを好みます。

これはあなたの「ループなし」要件を満たしていませんが、私ならこうするでしょう:

for m = 1:C

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

end

この結果、 A×D×C マトリックス Z.

そしてもちろん、次のように Z を事前に割り当てて処理を高速化することもできます。 Z = zeros(A,D,C);.

これは 1 行の解決策です (3 次元に分割する場合は 2 行)。

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

説明：

上記は複雑に見えるかもしれませんが、考え方は簡単です。まず、次の 3 次元を取得することから始めます。 X そして最初のディムに沿って垂直連結を行います。

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

...難しかったのは C は変数であるため、を使用してその式を一般化することはできません猫または バートキャット. 。次にこれを掛けます Y:

ZZ = XX * Y;

最後に、それを 3 次元に分割して戻します。

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

したがって、必要な行列の乗算は 1 つだけであることがわかりますが、次のことを行う必要があります。 形を変える 前後のマトリックス。

私は最も効率的な方法を念頭に置いて、まったく同じ問題に取り組んでいます。外部ライブラリを使用する以外に、私がよく目にするアプローチは大まかに 3 つあります (つまり、 mtimesx):

3D マトリックスのスライスをループします。
repmat と 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

3 つのアプローチはすべて同じ出力を生成しましたが (ふぅ!)、驚くべきことに、ループが最も高速でした。

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) が大きな要因で勝ちます。すべての場合に最適なアプローチは存在しない可能性がありますが、ループは依然として非常に優れており、多くの場合に最適です。読みやすさの点でも最高です。ループしてください！

いいえ。方法はいくつかありますが、直接的または間接的に常にループして発生します。

私の好奇心を満足させるために、そもそもなぜそれが欲しいのですか？

質問に答えると、 そして 読みやすさについては、以下を参照してください。

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

を使用することを強くお勧めします MMX ツールボックスマットラボの。n 次元行列を可能な限り高速に乗算できます。

の利点 MMX は：

それは簡単使用します。
かける n次元行列 (実際には 2 次元行列の配列を乗算できます)
その他を実行します 行列演算 (転置、二次乗算、Chol 分解など)
それは使用しています 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  <===================

私は再帰を考えますが、それがあなたが実行できる他の非ループメソッドだけです

ループを「展開」することもできます。つまり、ループ内で発生するすべての乗算を順番に書き出すことができます。

次の問題に対する私の答えを共有したいと思います。

1) 2 つのテンソル (任意の価数) のテンソル積を作成します。

2) 任意の次元に沿って 2 つのテンソルを縮小します。

最初と 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

批評家はいますか？

ライセンス： CC-BY-SA と帰属

所属していません StackOverflow