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

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

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

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

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

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

**追記（2020/09/18)**
- R2020b からは pagemtimes 関数（[ヘルプページ](https://www.mathworks.com/help/matlab/ref/pagemtimes.html?s_eid=PSM_29435)）が使えます。



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

```matlab
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: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 ループがあまり美しくない。使わずキレイなコードを書けないのか・・・

もちろん

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

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


# 試したこと＋その結果

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

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

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

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

```matlab: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 関数の生成](https://jp.mathworks.com/help/coder/ug/how-to-generate-mex-functions-using-the-matlab-coder-project-interface.html?s_eid=PSM_29435)

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

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

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

[File Exchange](https://jp.mathworks.com/matlabcentral/fileexchange?s_eid=PSM_29435) ってご存知ですか？ ユーザーが作った関数や、MATLAB 製品に入っていない機能が沢山公開されているページです。誰でも公開可能なのでよかったら作品を公開してみてください。

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

1. [MTIMESX - Fast Matrix Multiply with Multi-Dimensional Support](https://jp.mathworks.com/matlabcentral/fileexchange/25977-mtimesx-fast-matrix-multiply-with-multi-dimensional-support?s_eid=PSM_29435)[^3]
2. [MMX: N-D multithreaded matrix operations]
(https://jp.mathworks.com/matlabcentral/fileexchange/37515-mmx?s_eid=PSM_29435)[^4]

[^3]: 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.
[^4]: Yuval (2020). mmx (https://www.github.com/yuvaltassa/mmx), GitHub. Retrieved January 7, 2020.


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

[^1]: 詳細はこちら：[C MEX ファイル アプリケーション](https://jp.mathworks.com/help/matlab/call-mex-files-1.html?s_eid=PSM_29435)
[^2]: MinGW-64 はこちらから：[MATLAB Support for MinGW-w64 C/C++ Compiler](https://jp.mathworks.com/matlabcentral/fileexchange/52848-matlab-support-for-mingw-w64-c-c-compiler?s_eid=PSM_29435)


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

```matlab

mex -setup
```

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

### [MTIMESX - Fast Matrix Multiply with Multi-Dimensional Support](https://jp.mathworks.com/matlabcentral/fileexchange/25977-mtimesx-fast-matrix-multiply-with-multi-dimensional-support?s_eid=PSM_29435)

![mtimesx.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/149511/b36ae02a-706b-bb75-19cc-a98f61242909.png)

> 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用）

```matlab
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` ができますので、


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

てな感じで使えます。

### MMX: N-D multithreaded matrix operations

![mmx.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/149511/5a3ff557-2943-6ee1-7caa-9656ed3786e1.png)
> 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`

ができます。

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

と実行。

# まとめ

いろいろと試してみましたが、結局は [File Exchange](https://jp.mathworks.com/matlabcentral/fileexchange?s_eid=PSM_29435) に助けられました。さすが 36,000 を超える数のファイルがアップロードされているだけあって、何らかのものが見つかります。

今回登場した mex 関数については以下の記事も参考になります。
- [Qiita: mex function ことはじめ　～N次元配列となかよし編～](https://qiita.com/iyselee07/items/234e7dd362031222f5bf) by @iyselee07 さん
- [Qiita: [MATLAB] MEX関数の作成方法](https://qiita.com/shohirose/items/0689a76638fa217c6e83) by @shohirose さん

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

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