行列の積でベンチマーク : NumPy, Numba, Cython, Swig, OpenCL, intelMKL

\[2020-08-26追記\]
広告になりますが記事が古くなってきたのでコレスキー分解をベンチの対象として更新しました。
[Pythonと関数単位での高速化\(cython, numba, swig\) \- 夜間飛行](https://pid123.blogspot.com/2020/06/python-cython-numba-swig.html)
\[2020-08-26追記ここまで\]

組み込みクラスの環境にディープラーニングを突っ込むことに興味があるのですが、その中で実装も簡単で実行速度もよい方法を模索しています。
そこで、ちょっとしたベンチマークのため、NumPy, Numba, Cython, Swig, OpenCL, intelMKLでの **NxN** 行列2つの積を計算させました。

Swigは中でfloatをdoubleに非明示的にキャストしてたり、実装確認中のメモが一部ゾンビ化してしまいましたが、ご査収ください。もっと比較水準増やしたり、環境違いの結果なPL待ってまーす。

[NumPyCuPy\_test/Test01\.ipynb at master · Chachay/NumPyCuPy\_test](https://github.com/Chachay/NumPyCuPy_test/blob/master/Test01.ipynb)

# 今回掛け算される人

単精度のNxN行列 maと同じく単精度のNxN行列 mb

```python:行列の定義
N = 1000 # 下の表のNと一致
# Generate NxN randomized matrix
ma = np.random.rand(N, N).astype(np.float32)
mb = np.random.rand(N, N).astype(np.float32)
```

# 結果まとめ

NumbaでNumpyのdotを使った結果が圧倒的。だけど、実行結果キャッシュされてたりしない？

|   |                              | N=3          |      | N=100        |      | N=1000       |      |
|---|------------------------------|--------------|------|--------------|------|--------------|------|
|   |                              | 実行時間[us] | 速度 | 実行時間[us] | 速度 | 実行時間[ms] | 速度 |
| 1   | Numpy Native                 |         1.34 | 1.00 |         94.2 | 1.00 |         87.9 | 1.00 |
| 2   | Numba Simple Mult            |         1.24 | 1.08 |      1,040.0 | 0.09 |      6,870.0 | 0.01 |
| 3   | Numba Numpy                  |         2.03 | 0.66 |         74.6 | 1.26 |         23.5 | 3.74 |
| 4   | Cython Simple Mult           |        19.20 | 0.07 |    660,000.0 | 0.00 |           NA |    - |
| 5   | Cython Numpy                 |         1.52 | 0.88 |         94.1 | 1.00 |         87.8 | 1.00 |
| 6   | Swig Simple Mult             |         1.97 | 0.68 |        837.0 | 0.11 |      7,450.0 | 0.01 |
| 7   | Swig MKL cblas_dgemm         |         2.57 | 0.52 |        123.0 | 0.77 |         95.8 | 0.92 |
| 8   | OpenCL Simple Mult(Parallel) |     2,560.00 | 0.00 |      3,160.0 | 0.03 |        877.0 | 0.10 |

[メモ] 実行しながらメモったはずだけど大小関係おかしい。みなおしてみます。

# Numpy

リファレンス見てください

```python
np.dot(ma, mb)
```

# Numba

np.dotにjitをかけたNum_NpDotと、数学の定義通り行列掛け算をするNum_Dotを準備

```python
from numba import jit
@jit
def Num_NpDot(a, b):
    return np.dot(a, b)
@jit    
def Num_Dot(a, b):
    c = np.zeros((a.shape[0], b.shape[1]))
    for i in range(a.shape[0]):
        for j in range(b.shape[1]):
            for k in range(a.shape[1]):
                c[i][j] += a[i][k] * b[k][j]
    return c
```

# Cython

型の定義もしてやるなど目をかけながら世話をして、Cに翻訳しているということからNumbaには勝てるだろうと思っていたら、今回のような結果になり、お前にはかなりがっかりしたぞ。

Numba水準と同じ. ビルドはVC2015でO2相当の最適化（のはず).

```python
cimport numpy as np
cimport cython
import numpy as np

@cython.boundscheck(False) # turn off bounds-checking for entire function
@cython.wraparound(False)  # turn off negative index wrapping for entire function

cpdef np.ndarray Cy_NpDot(np.ndarray a, np.ndarray b):
    return np.dot(a, b)
    
cpdef np.ndarray Cy_Dot(np.ndarray a, np.ndarray b):
    cdef np.ndarray c
    c = np.zeros((a.shape[0], b.shape[0]))
    for i in range(a.shape[0]):
        for j in range(b.shape[1]):
            for k in range(a.shape[1]):
                c[i][j] += a[i][k] * b[k][j]
    return c
```

Cy_DotがCython通しているくせにとても使えない子になって困った。

# Swig

最近すっかりお世話になるようになった。PythonやってるのにC++で大事なところを書く。
Pythonはデータハンドリングにしか使わないのですよ、方式。

Numbaと同じく、いわゆる掛け算のSwig_DotとNumpy相当水準としてIntelMKL登板。
ビルドはVC2015でO2相当の最適化（のはず).

単精度なのにdoubleで用意してしまった。しまった。
(特にMKLの中で変なことになっていそうだけど、とりあえず目をつぶる)

```C++
void Swig_Dot(int mm1, int mn1, double* mat1,
              int mm2, int mn2, double* mat2,
              double** outmat, int* mm, int* mn)
{
    double* arr = NULL;
    arr = (double*)calloc(mm1 * mn2, sizeof(double));
    
    // Multiplying matrix a and b and storing in array mult.
    for(int i = 0; i < mm1; ++i)
        for(int j = 0; j < mn2; ++j)
            for(int k = 0; k < mn1; ++k)
                arr[i*mn2+j] += mat1[i*mn1+k] * mat2[k*mn2+j];
    *mm = mm1;
    *mn = mn2;
    *outmat = arr;
}

void Swig_Dot_MKL(int mm1, int mn1, double* mat1,
              int mm2, int mn2, double* mat2,
              double** outmat, int* mm, int* mn)
{
    double alpha = 1.0;
    double beta = 0.0;

    double *C = (double *)mkl_malloc( mm1*mn2*sizeof( double ), 64 );
    if (C == NULL) {
        printf( "\n ERROR: Can't allocate memory for matrices. Aborting... \n\n");
        mkl_free(C);
        return;
    }

    for (int i = 0; i < (mm1*mn2); i++) {
        C[i] = 0.0;
    }
    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, 
                mm1, mn2, mn1, alpha, mat1, mn1, mat2, mn2, beta, C, mn2);

    double *res = (double *)malloc( mm1*mn2*sizeof( double ));
    memcpy(res, C, mm1*mn2*sizeof( double ));
    *mm = mm1;
    *mn = mn2;
    *outmat = res;
    mkl_free(C);
}
```

# OpenCL

おまけ。GPUと比較するためのスタブだと思ってください。CuPyと比べようと思ったけど、ただのハードウェア マッチョクラブになりそう。

省略。


# 参考文献
* [Python \(NumPy\) と Common Lisp \(LLA\) で行列積の実行速度を比較する \- 不確定特異点](http://tanakahx.hatenablog.com/entry/2015/09/25/070000)
* [Parallel programming with opencl and python – prototype](http://www.nehalemlabs.net/prototype/blog/2014/04/28/parallel-programming-with-opencl-and-python/)
* [10\.3: OpenCL](http://bathompso.com/education/compphys/parallel/opencl/)

* [Multiplying Matrices Using dgemm \| Intel® Software](https://software.intel.com/en-us/node/529735)
* [Intel® Math Kernel Library Link Line Advisor \| Intel® Software](https://software.intel.com/en-us/articles/intel-mkl-link-line-advisor)
* [cblas\_dgemm \- Accelerate \| Apple Developer Documentation](https://developer.apple.com/reference/accelerate/1513282-cblas_dgemm?language=objc)
* [PyOpenCLによるGPGPU入門](https://www.slideshare.net/likr/pyopenclgpgpu)
* [どうすればPythonをJuliaと同じくらい速く動かせるのか？ : 様々なやり方で計算の高速化を図る \| プログラミング \| POSTD](http://postd.cc/python_meets_julia_micro_performance/)