[2020-08-26追記]
広告になりますが記事が古くなってきたのでコレスキー分解をベンチの対象として更新しました。
Pythonと関数単位での高速化(cython, numba, swig) - 夜間飛行
[2020-08-26追記ここまで]
組み込みクラスの環境にディープラーニングを突っ込むことに興味があるのですが、その中で実装も簡単で実行速度もよい方法を模索しています。
そこで、ちょっとしたベンチマークのため、NumPy, Numba, Cython, Swig, OpenCL, intelMKLでの NxN 行列2つの積を計算させました。
Swigは中でfloatをdoubleに非明示的にキャストしてたり、実装確認中のメモが一部ゾンビ化してしまいましたが、ご査収ください。もっと比較水準増やしたり、環境違いの結果なPL待ってまーす。
NumPyCuPy_test/Test01.ipynb at master · Chachay/NumPyCuPy_test
今回掛け算される人
単精度のNxN行列 maと同じく単精度のNxN行列 mb
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
リファレンス見てください
np.dot(ma, mb)
Numba
np.dotにjitをかけたNum_NpDotと、数学の定義通り行列掛け算をするNum_Dotを準備
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相当の最適化（のはず).
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の中で変なことになっていそうだけど、とりあえず目をつぶる)
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と比べようと思ったけど、ただのハードウェア マッチョクラブになりそう。
省略。