3

forループを排除して速度を上げるために、「ベクトル化」された方法で、大きな（2x2m）マトリックスの（2x2）サブマトリックスを乗算しようとしています。現在、(2x2xm) に変形してから、for ループを使用してこれを行います。

for n = 1:1e5
    m = 1e4;
    A = rand([2,2*m]);     % A is a function of n
    A = reshape(A,2,2,[]);
    B = eye(2);
    for i = 1:m
        B = A(:,:,i)*B;    % multiply the long chain of 2x2's
    end
end

関数のゴールは @prod に似ていますが、要素ごとのスカラー乗算ではなく行列乗算を使用します。@multiprod は近いように見えますが、引数として 2 つの異なる nD 行列を取ります。非常に大きな 2D 配列の複数の部分行列、または単一の 2x2m{xn} 配列を使用して for ループの 1 つまたは両方を排除するソリューションを想像します。

前もって感謝します、ジョー

4

2 に答える 2

0

以下のコードのように、ベクトル化された乗算を行うには、別の方法で行列を再形成する必要があると思います。このコードもループを使用していますが、より高速である必要があると思います

MM      = magic(2);
M0      = MM;
M1      = rot90(MM,1);
M2      = rot90(MM,2);
M3      = rot90(MM,3);


MBig1           = cat(2,M0,M1,M2,M3);
fprintf('Original matrix\n')
disp(MBig1)
MBig2           = zeros(size(MBig1,2));
MBig2(1:2,:)    = MBig1;
for k=0:3
    c1 =  k   *2+1;
    c2 = (k+1)*2+0;
    MBig2(:,c1:c2) = circshift(MBig2(:,c1:c2),[2*k 0]);
end
fprintf('Reshaped original matrix\n')
disp(MBig2)

fprintf('Checking [ M0*M0 M0*M1 M0*M2 M0*M3 ] in direct way\n')
disp([ M0*M0 M0*M1 M0*M2 M0*M3 ])
fprintf('Checking [ M0*M0 M0*M1 M0*M2 M0*M3 ] in vectorized way\n')
disp( kron(eye(4),M0)*MBig2 )


fprintf('Checking [ M0*M1*M2*M3 ] in direct way\n')
disp([ M0*M1*M2*M3 ])
fprintf('Checking [ M0*M1*M2*M3 ] in vectorized way\n')
R2 = MBig2;
for k=1:3
    R2 = R2 * circshift(MBig2,-[2 2]*k);
end
disp(R2)

出力は

Original matrix
     1     3     3     2     2     4     4     1
     4     2     1     4     3     1     2     3

Reshaped original matrix
     1     3     0     0     0     0     0     0
     4     2     0     0     0     0     0     0
     0     0     3     2     0     0     0     0
     0     0     1     4     0     0     0     0
     0     0     0     0     2     4     0     0
     0     0     0     0     3     1     0     0
     0     0     0     0     0     0     4     1
     0     0     0     0     0     0     2     3

Checking [ M0*M0 M0*M1 M0*M2 M0*M3 ] in direct way
    13     9     6    14    11     7    10    10
    12    16    14    16    14    18    20    10

Checking [ M0*M0 M0*M1 M0*M2 M0*M3 ] in vectorized way
    13     9     0     0     0     0     0     0
    12    16     0     0     0     0     0     0
     0     0     6    14     0     0     0     0
     0     0    14    16     0     0     0     0
     0     0     0     0    11     7     0     0
     0     0     0     0    14    18     0     0
     0     0     0     0     0     0    10    10
     0     0     0     0     0     0    20    10

Checking [ M0*M1*M2*M3 ] in direct way
   292   168
   448   292

Checking [ M0*M1*M2*M3 ] in vectorized way
   292   168     0     0     0     0     0     0
   448   292     0     0     0     0     0     0
     0     0   292   336     0     0     0     0
     0     0   224   292     0     0     0     0
     0     0     0     0   292   448     0     0
     0     0     0     0   168   292     0     0
     0     0     0     0     0     0   292   224
     0     0     0     0     0     0   336   292
于 2013-06-21T06:59:47.047 に答える
0

以下の関数は、私の問題の一部を解決するかもしれません。これは、times と mtimes のように、"mprod" と prod という名前が付けられています。いくつかの再形成により、multiprodを再帰的に使用します。一般に、再帰関数呼び出しはループよりも低速です。Multiprod は 100 倍以上高速であると主張しているため、それを補って余りあるはずです。

function sqMat = mprod(M)
    % Multiply *many* square matrices together, stored
    % as 3D array M. Speed gain through recursive use 
    % of function 'multiprod' (Leva, 2010).

    % check if M consists of multiple matrices
    if size(M,3) > 1
        % check for odd number of matrices
        if mod(size(M,3),2)
            siz = size(M,1);
            M = cat(3,M,eye(siz));
        end
        % create two smaller 3D arrays
        X = M(:,:,1:2:end); % odd pages
        Y = M(:,:,2:2:end); % even pages
        % recursive call
        sqMat = mprod(multiprod(X,Y));
    else
        % create final 2D matrix and break recursion
        sqMat = M(:,:,1);
    end
end

この機能の速度や精度はテストしていません。 これはループよりもはるかに高速だと思います。高次元では使用できないため、操作を「ベクトル化」しません。この関数を繰り返し使用する場合は、ループ内で行う必要があります。

EDIT以下は、十分に高速に動作するように見える新しいコードです。関数の再帰呼び出しは遅く、スタック メモリを消費します。まだループを含んでいますが、log(n)/log(2) だけループの数を減らします。また、より多くのディメンションのサポートが追加されました。

function sqMats = mprod(M)
    % Multiply *many* square matrices together, stored along 3rd axis.
    % Extra dimensions are conserved; use 'permute' to change axes of "M".
    % Speed gained by recursive use of 'multiprod' (Leva, 2010).

    % save extra dimensions, then reshape
    dims = size(M);
    M = reshape(M,dims(1),dims(2),dims(3),[]);
    extraDim = size(M,4);

    % Check if M consists of multiple matrices...
    % split into two sets and multiply using multiprod, recursively
    siz = size(M,3);
    while siz > 1
        % check for odd number of matrices
        if mod(siz,2)
            addOn = repmat(eye(size(M,1)),[1,1,1,extraDim]);
            M = cat(3,M,addOn);
        end
        % create two smaller 3D arrays
        X = M(:,:,1:2:end,:); % odd pages
        Y = M(:,:,2:2:end,:); % even pages
        % recursive call and actual matrix multiplication
        M = multiprod(X,Y);
        siz = size(M,3);
    end

    % reshape to original dimensions, minus the third axis.
    dims(3) = [];
    sqMats = reshape(M,dims);
end
于 2013-06-22T03:04:42.460 に答える