質問: 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文で回せ」という回答は最もだと思いますが、理解しやすさはここでは不問とさせて下さい。
回答: 下記が参考になるかもしれません。 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 ご参考まで。