N-D array（N次元配列）の計算

先週とみー @tommyecguitar さんからこんな質問頂きました。

多次元配列の1、2次元を行列とみなした乗算ってできないんですかね？for 使うしかないんですかね。
(Qiita: このプロットどうやって描いたの？のコメント欄参照。)

例えば 20x20x1000 と 20x20x1000 の配列があったときに 20x20 の行列演算を 1000 回実行する場合です

結論から言うと

MATLAB 関数使う限りは for ループが最も効率よさそう
C-mex 関数を使うと for ループの 4 倍程度高速に計算できました

追記（2020/09/18)

R2020b からは pagemtimes 関数（ヘルプページ）が使えます。

コードとそれぞれの実行時間は以下の通り。

a = rand(20,20,1000);
b = rand(20,20,1000);

% for ループでの計算（0.0057 秒）
tic, c = matlab_mprod(a,b); toc 

% mex 関数化して計算（0.0046 秒）
tic, c = matlab_mprod_mex(a,b); toc

% mtimesx 関数利用（0.0026 秒）
tic, c = mtimesx(a,b); toc

% mmx 関数利用（0.0017 秒）
tic, c = mmx('mult',a,b); toc

といった具合。R2019b (Windows 10) で試しました。

まず for ループではイケナイの？

まず思いつく方法としては for ループを使った

matlab_mprod.m


function c = matlab_mprod(a,b)
c  = zeros(size(a,1),size(b,2),size(a,3));
for i = 1:size(a,3)
   c(:,:,i) = a(:,:,i)*b(:,:,i);
end

ですが、for ループがあまり美しくない。使わずキレイなコードを書けないのか・・・

もちろん

c = cell2mat(cellfun(@(x,y) x*y, num2cell(a,[1 2]), num2cell(b,[1 2]), 'UniformOutput', false));

のように cellfun や num2cell を使って 1 行で書くことは可能なんですが、これは計算速度が遅い。for ループ使った方法の 4 倍程度時間がかかります。

試したこと＋その結果

スパース行列の利用（遅い）
MATLAB Coder による mex 関数化（多少速い）
File Exchange で公開されている関数の使用（速い）

スパース行列の利用（遅い）

狙い：20x20x1000 の 2 次元配列を 2 次元の配列にしてやって、掛算を 1 回だけの実行にすれば速いのでは？
結果：めっちゃ遅い

200000x200000 のブロック対角行列（スパース）の掛け算を実施することになるのですが、スパース行列の演算を利用すれば速くなるのではという期待がありましたが、見事に裏切られました。

ndprod_sparse.m

N = 1000;
a = rand(20,20,N);
b = rand(20,20,N);
aa = num2cell(a,[1,2]);
bb = num2cell(b,[1,2]);
aa{1} = sparse(aa{1});
bb{1} = sparse(bb{1});
A = blkdiag(aa{:});
B = blkdiag(bb{:});
C = A*B;

てな感じです。計算結果 C を 3 次元配列に直す処理も考えると絶望的ですね。for ループ使った方がコードも短いので、もはや救いようがありません。すいません。

MATLAB Coder による mex 関数化（多少速い）

狙い：C コード生成しちゃえば速いのでは？
結果：MATLAB の for ループより多少速い程度（R2019b）

R2015b から MATLAB の実行エンジンが刷新されて、JIT (Just in Time) コンパイルにより for ループの実行自体の処理速度はかなり改善しました。これはむしろ MATLAB の for ループすごいじゃん・・とするべきか。

ちなみに mex 関数というのは、C や Fortran で書かれたコードから MATLAB から実行できる形にコンパイルされた関数のことを言います。MATLAB Coder を使うと m -> C/C++ -> mex 関数の変換を自動化できる。詳細は MATLAB Coder アプリを使用した MEX 関数の生成

ただ、多少速くなったとはいえ安くはない MATLAB Coder を使ってこの結果は納得がいきません。

File Exchange で公開されている関数の使用（速い）

狙い：もう既に誰かが速い関数作っているんじゃない？
結果：すごいのあった

File Exchange ってご存知ですか？ユーザーが作った関数や、MATLAB 製品に入っていない機能が沢山公開されているページです。誰でも公開可能なのでよかったら作品を公開してみてください。

ここで見つけたのが以下の 2 つ。

MTIMESX - Fast Matrix Multiply with Multi-Dimensional Support¹
[MMX: N-D multithreaded matrix operations]
(https://jp.mathworks.com/matlabcentral/fileexchange/37515-mmx?s_eid=PSM_29435)[^4]

今回は使ってみただけですが、どちらもかなり気合の入ったもので、今回の使い方以外にも沢山機能があります。興味のある方は深堀りしてみてください。

どちらも C コードをコンパイルして mex 関数²を作る必要があるので C のコンパイラが必要です。無料の MinGW-64³でも大丈夫です。インストール後に


mex -setup

でコンパイラを設定してください。使うまでのポイントだけ紹介します。

MTIMESX - Fast Matrix Multiply with Multi-Dimensional Support

MTIMESX is a fast general purpose matrix and scalar multiply routine that has the following features:

Supports multi-dimensional (nD, n>2) arrays directly
Supports Transpose, Conjugate Transpose, and Conjugate pre-operations
Supports singleton expansion
Utilizes BLAS calls, custom C loop code, or OpenMP multi-threaded C loop code
Can match MATLAB results exactly or approximately as desired
Can meet or beat MATLAB for speed in most cases

とのこと。紹介文だけで速そうな感じしますね。

リンク先からファイルをダウンロードすると、mtimesx.c があります。これをコンパイルして mtimesx.mexw64 というバイナリを作成します。コンパイルするための関数 mtimesx_build.m も提供されているのですが、私の環境ではうまく動きませんでした。

ですので、File Exchange のページのコメントにある以下のコードでコンパイルします。（Windows用）

libdir = 'microsoft';
comp = computer;
mext = mexext;
lc = length(comp);
lm = length(mext);
cbits = comp(max(1:lc-1):lc);
mbits = mext(max(1:lm-1):lm);
if( isequal(cbits,'64') || isequal(mbits,'64') )
compdir = 'win64';
largearraydims = '-largeArrayDims';
else
compdir = 'win32';
largearraydims = '';
end
lib_blas = [matlabroot '\extern\lib\' compdir '\' libdir '\libmwblas.lib'];
d = dir('mtimesx.m');
mname = [d.folder '/' d.name];
cname = [mname(1:end-2) '.c'];
mex(cname,largearraydims,lib_blas,'COMPFLAGS="$COMPFLAGS /openmp"');

これを実行するとmtimesx.mexw64 ができますので、

a = rand(20,20,1000);
b = rand(20,20,1000);
c = mtimesx(a,b);

てな感じで使えます。

MMX: N-D multithreaded matrix operations

Multithreaded matrix operations on N-D arrays (a Matlab plug-in)

だそうな。こちらはファイルをダウンロード後に src フォルダ内で build_mmx.m を実行するだけでOK.

mmx.mexw64
mmx_mkl_multi.mexw64
mmx_naive.mex64

ができます。

a = rand(20,20,1000);
b = rand(20,20,1000);
c = mmx(a,b);

と実行。

まとめ

いろいろと試してみましたが、結局は File Exchange に助けられました。さすが 36,000 を超える数のファイルがアップロードされているだけあって、何らかのものが見つかります。

今回登場した mex 関数については以下の記事も参考になります。

Qiita: mex function ことはじめ　～N次元配列となかよし編～ by @iyselee07 さん
Qiita: [MATLAB] MEX関数の作成方法 by @shohirose さん

今回調査する機会を提供してくださったとみー@tommyecguitar さん、ありがとうございました！

他の方も何か気になることがあればコメントお待ちしています。

James Tursa (2020). MTIMESX - Fast Matrix Multiply with Multi-Dimensional Support (https://www.mathworks.com/matlabcentral/fileexchange/25977-mtimesx-fast-matrix-multiply-with-multi-dimensional-support?s_eid=PSM_29435), MATLAB Central File Exchange. Retrieved January 7, 2020. ↩
詳細はこちら：C MEX ファイルアプリケーション ↩
MinGW-64 はこちらから：MATLAB Support for MinGW-w64 C/C++ Compiler ↩

【MATLAB】N-D array（N次元配列）の計算速度比較