2016-07-28 8 views
3

ここでは、元のコードです:ベクトル化およびネストされた行列乗算

K = zeros(N*N) 
for a=1:N 
    for i=1:I 
     for j=1:J 
      M = kron(X(:,:,a).',Y(:,:,a,i,j)); 

      %A function that essentially adds M to K. 
     end 
    end 
end 

目標はkroniker乗算呼び出しをベクトル化することです。私の直感は、XとYを行列のコンテナと考えることです(参考のために、kronに供給されるXとYのスライスは7x7の正方行列です)。このコンテナスキームでは、Xは1次元のコンテナとして表示され、Yは3次元のコンテナとして表示されます。私の次の推測は、Yを2次元のコンテナに変換するか、または1次元のコンテナに変換し、次にXとYを要素的に掛け合わせることでした。質問:どのようにしてこの変形をMとMATLABはこのコンセプトでこのアイデアを扱うことができますか?また、内側の行列要素をさらに公開するためにコンテナをさらに再構成する必要がありますか?

答えて

6

アプローチ#1:シンプル7D並び替える

% Get sizes 
[m1,m2,~] = size(X); 
[n1,n2,N,n4,n5] = size(Y); 

% Perform kron format elementwise multiplication betwen the first two dims 
% of X and Y, keeping the third dim aligned and "pushing out" leftover dims 
% from Y to the back 
mults = bsxfun(@times,permute(X,[4,2,5,1,3]),permute(Y,[1,6,2,7,3,4,5])); 

% Lose the two dims with summation reduction for final output 
out = sum(reshape(mults,m1*n1,m2*n2,[]),3); 

検証

6D並び替える

% Get sizes 
[m1,m2,~] = size(X); 
[n1,n2,N,n4,n5] = size(Y); 

% Lose the third dim from X and Y with matrix-multiplication 
parte1 = reshape(permute(Y,[1,2,4,5,3]),[],N)*reshape(X,[],N).'; 

% Rearrange the leftover dims to bring kron format 
parte2 = reshape(parte1,[n1,n2,I,J,m1,m2]); 

% Lose dims correspinding to last two dims coming in from Y corresponding 
% to the iterative summation as suggested in the question 
out = reshape(permute(sum(sum(parte2,3),4),[1,6,2,5,3,4]),m1*n1,m2*n2) 

アプローチ#2と行列の乗算我々は最高を参照してください、実行した後

% Setup inputs 
X = rand(10,10,10); 
Y = rand(10,10,10,10,10); 

% Original approach 
[n1,n2,N,I,J] = size(Y); 
K = zeros(100); 
for a=1:N 
    for i=1:I 
     for j=1:J 
      M = kron(X(:,:,a).',Y(:,:,a,i,j)); 
      K = K + M; 
     end 
    end 
end 

% Approach #1 
[m1,m2,~] = size(X); 
[n1,n2,N,n4,n5] = size(Y); 
mults = bsxfun(@times,permute(X,[4,2,5,1,3]),permute(Y,[1,6,2,7,3,4,5])); 
out1 = sum(reshape(mults,m1*n1,m2*n2,[]),3); 

% Approach #2 
[m1,m2,~] = size(X); 
[n1,n2,N,n4,n5] = size(Y); 
parte1 = reshape(permute(Y,[1,2,4,5,3]),[],N)*reshape(X,[],N).'; 
parte2 = reshape(parte1,[n1,n2,I,J,m1,m2]); 
out2 = reshape(permute(sum(sum(parte2,3),4),[1,6,2,5,3,4]),m1*n1,m2*n2); 

から

はここで、オリジナルと提案したアプローチを実行するための設定です。元のものに対する提案されたアプローチとの絶対偏差 -

>> error_app1 = max(abs(K(:)-out1(:))) 
error_app1 = 
    1.1369e-12 
>> error_app2 = max(abs(K(:)-out2(:))) 
error_app2 = 
    1.1937e-12 

値は私によく見えます!


ベンチマーク

検証に使用され、私たちはこのような何かを得るのと同じ大きなデータセットを使用して、これらの3つのアプローチのタイミング -

----------------------------- With Loop 
Elapsed time is 1.541443 seconds. 
----------------------------- With BSXFUN 
Elapsed time is 1.283935 seconds. 
----------------------------- With MATRIX-MULTIPLICATION 
Elapsed time is 0.164312 seconds. 

をマトリックス乗算はかなり良いやっているように思えますこれらのサイズのデータ​​セットに対して!

+0

新しいアプローチの両方で、out1とout2はKを置き換えますか? もしそうなら、K = K + Mがより複雑な場合、これを反映することは可能でしょうか? たとえば、元の例では、各Mがループの反復ごとに変化する定数で正規化された場合はどうなりますか? 2番目の例は、Kに追加する前にMにブール値マスクが適用されている場合です。これらの余分なステップで出力が中断されますか? –

+0

@MikeVandenberg前処理ステップとして、 'Y = bsxfun(@ times、Y、permute(scale、[4,5,1,2,3]));を使用し、提案コードを使用します。 scaleはサイズ '(N、I、J)'のスケーリング行列です。 ex2の場合は、すべての反復または同じ反復に対して異なるマスクを使用しますか? – Divakar

+0

私はこれらのステップについて単純化してはいけないと思います。以下は、Kを構成する特定の関数です: 'pr = real(trace(E * M)); ここで、Hは、手前の既知の行列です。各反復ではいですので、Hの別のスライスを引っ張って、Mを正規化すると、積E * Mのトレースである異なる定数で正規化します。ここで、Eはブール値のマスクです。 –

関連する問題