MATLAB Answers

Daichi
0

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

Asked by Daichi
on 22 Mar 2019
Latest activity Edited by Tohru Kikawada on 28 Apr 2019
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 Comment

SIMD命令と言われている内容ですが、特定のプロセッサにおけるベクトル演算ユニット等への命令と理解しましたが、よろしいでしょうか。仮にそうだとして、「コロン演算子を用いる書き方ではSIMD命令が使われるから」というのが正しいのかは良くわからないです 。
なぜなら、「for文でNに関してループしつつdetやinvするよりはるかに高速です。」とありますが、MATLABのdetやinvは、
にあるアルゴリズム(LU分解)を使っており内部で演算精度を保持するため、ピボット選択を行なっていますから、条件判断分岐が入ります。このような演算を行わずに公式をそのまま使った計算では早くなるのが当然と思われます。Mが大きくなると更に数値計算の誤差を考慮しない計算は危険だと思います。またこれは想像ですが、detやinv自体の計算では、恐らくベクトル演算ユニットが使えるようにCPUに応じた命令セットを使っているのではと考えます。

Sign in to comment.

Products


Release

R2019a

1 Answer

Answer by Tohru Kikawada on 28 Apr 2019
Edited by Tohru Kikawada on 28 Apr 2019

下記が参考になるかもしれません。
基本的には任意のM対して下記を参考にご自身でベクトル化を施す必要があります。
ご指摘のとおり、項数が多くなるような場合には手計算や手書きでコードを書くのは大変です。下記のように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(:))))
計算結果:
ループで計算した場合の処理時間
経過時間は 0.077265 秒です。
ベクトル化した場合の処理時間
経過時間は 0.002525 秒です。
要素ごとの差分の最大値:
2.2880e-07
詳細は下記をご覧ください。
また、パフォーマンス面に関して、MATLABは線形代数計算は実行されているシステム(OS, CPU)のアーキテクチャに最適化されたライブラリを使って解いています。
その中ではSIMD演算などのCPUの並列専用命令の利用やメモリ帯域の最適化など様々な工夫が行われています。
このケースでは計算精度は考慮しないということかと思いますが、Matrix Vector様のコメントにもあるように逆行列を直接求めることは計算誤差の影響が無視できなくなるケースもありますので、注意が必要です。
ご参考まで。

  0 Comments

Sign in to comment.