2次元行列演算のSIMD命令による高速化

Question

Daichi 2019 年 3 月 22 日

0
リンク

この質問への直接リンク

https://jp.mathworks.com/matlabcentral/answers/451698-2-simd

編集済み: Tohru Kikawada 2019 年 4 月 28 日

M=2;

N=10000;

A=rand(M,M,N);

のようなサイズ感のAに対して、MxMの2次元行列の行列式や逆行列をN個全てについて求めるのはどうすれば高速でしょうか。 M=2の時はAの3次元目にコロン演算子を用いて

detA=A(1,1,:)*A(2,2,:)-A(1,2,:)*A(2,1,:);

などとすれば、for文でNに関してループしつつdetやinvするよりはるかに高速です。これは、上記のコロン演算子を用いる書き方ではSIMD命令が使われるからだと理解しています。 M=3の場合もサラスの法則を使えば同様のことができますが、Mが大きくなるとこれらの演算方法では項数が爆発的に増えていき、対処しきれなくなります。一般のMについて、SIMDを使えるような工夫があれば教えて下さい。 pythonではアインシュタインの縮約記法のeinsumを使えば、同様の演算を高速にできるという話を聞いたことがあります（詳しくはありませんが）。MATLABで同様のことはできますか？

ただし、ここではとにかく高速化に重きをおいています。「理解し難いコードになるからfor文で回せ」という回答は最もだと思いますが、理解しやすさはここでは不問とさせて下さい。

1 件のコメント
-1 件の古いコメントを表示-1 件の古いコメントを非表示

Matrix Vector 2019 年 3 月 26 日

SIMD命令と言われている内容ですが、特定のプロセッサにおけるベクトル演算ユニット等への命令と理解しましたが、よろしいでしょうか。仮にそうだとして、「コロン演算子を用いる書き方ではSIMD命令が使われるから」というのが正しいのかは良くわからないです。

なぜなら、「for文でNに関してループしつつdetやinvするよりはるかに高速です。」とありますが、MATLABのdetやinvは、

https://jp.mathworks.com/help/matlab/ref/det.html

https://jp.mathworks.com/help/matlab/ref/inv.html

にあるアルゴリズム（LU分解）を使っており内部で演算精度を保持するため、ピボット選択を行なっていますから、条件判断分岐が入ります。このような演算を行わずに公式をそのまま使った計算では早くなるのが当然と思われます。Mが大きくなると更に数値計算の誤差を考慮しない計算は危険だと思います。またこれは想像ですが、detやinv自体の計算では、恐らくベクトル演算ユニットが使えるようにCPUに応じた命令セットを使っているのではと考えます。

サインインしてコメントする。

サインインしてこの質問に回答する。

Answer 1

Tohru Kikawada 2019 年 4 月 28 日

0
リンク

この回答への直接リンク

https://jp.mathworks.com/matlabcentral/answers/451698-2-simd#answer_372610

編集済み: Tohru Kikawada 2019 年 4 月 28 日

下記が参考になるかもしれません。

https://jp.mathworks.com/matlabcentral/answers/224361-simultaneously-inverting-many-matrices

基本的には任意のM対して下記を参考にご自身でベクトル化を施す必要があります。

https://jp.mathworks.com/help/matlab/matlab_prog/vectorization.html

ご指摘のとおり、項数が多くなるような場合には手計算や手書きでコードを書くのは大変です。下記のようにSymbolic Math Toolboxを使って解析解を求め、数値計算することもできます。

clear;
M=2;
N=10000;
A=rand(M,M,N);
% ループ計算の場合
invA_loop = zeros(size(A));
tic
for k = 1:N
    invA_loop(:,:,k) = inv(A(:,:,k));
end
disp('ループで計算した場合の処理時間');
toc
% ベクトル化した場合
As = sym('a', [M,M]); % シンボル式として定義
invAs = reshape(inv(As),[],1); % シンボルを逆行列で解く
invAfh = matlabFunction(invAs,'Vars',As); % 数値計算用の無名関数ハンドルに変換
tic
invA_sym = reshape(invAfh(A(1,1,:),A(2,1,:),A(1,2,:),A(2,2,:)),M,M,N);
disp('ベクトル化した場合の処理時間');
toc
% 要素ごとの差分の最大値
disp('要素ごとの差分の最大値：');
disp(max(abs(invA_loop(:)-invA_sym(:))))