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

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文で回せ」という回答は最もだと思いますが、理解しやすさはここでは不問とさせて下さい。

[Matrix Vector]

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]

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

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(:))))

計算結果：

ループで計算した場合の処理時間

経過時間は 0.077265 秒です。

ベクトル化した場合の処理時間

経過時間は 0.002525 秒です。

要素ごとの差分の最大値：

2.2880e-07

詳細は下記をご覧ください。

https://jp.mathworks.com/help/symbolic/inv.html

https://jp.mathworks.com/help/symbolic/matlabfunction.html

また、パフォーマンス面に関して、MATLABは線形代数計算は実行されているシステム(OS, CPU)のアーキテクチャに最適化されたライブラリを使って解いています。

https://jp.mathworks.com/help/matlab/math/lapack-in-matlab.html

その中ではSIMD演算などのCPUの並列専用命令の利用やメモリ帯域の最適化など様々な工夫が行われています。

このケースでは計算精度は考慮しないということかと思いますが、Matrix Vector様のコメントにもあるように逆行列を直接求めることは計算誤差の影響が無視できなくなるケースもありますので、注意が必要です。

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

ご参考まで。