CythonによるPythonの高速化を検証・分析してみた

# はじめに

動的型付けに分類される言語「Python」には、PythonコードをC/C++ に変換しコンパイルを行う「Cython」やJITコンパイラを使ってPythonを高速化するモジュールである「Numba」等があります。
これらを研究やインターンで利用しているPythonコードに上手く組み込むことで効率化を図りたいと考え、理解を深めるために実際の簡単な処理に対する処理速度の比較検証をjupyter上で行ってみました。

本記事はCythonの入門記事ではありません。Cythonとは何か？基礎的な知識、利用法等を知りたい場合は以下の記事をお薦めします。

- [cython入門](https://qiita.com/en3/items/1f1a609c4d7c8f3066a7)
- [深入りしないCython入門](https://qiita.com/pashango2/items/45cb85390193d97523ca)



意見や誤りがありましたらご気軽にコメントください。
実験コードは[コチラ](https://colab.research.google.com/drive/1Jeahw-_cdj4OWA2yfOtRLLq4o1t_u5qI)



# 実験1 is_prime関数での比較

以下に示される単純な`is_prime`関数（素数判定関数)をベースラインとして速度の高速化を試みました。コードについてベースラインと同じ箇所は略してあります。Cythonに関しては型付けの場所による速度の違いを検証しました。


##### baseline

```python

def find_factor(n):
    answer = n
    for i in range (2,n):
        if n % i == 0:
            answer = i
            break
    return answer

def is_prime(n):
    if find_factor (n) != n:
        return False
    else:
        return True
```

##### Numbaでの高速化
```python
from numba import jit
@jit
def find_factor0 (n):
    #略
def is_prime0 (n):
    #略
```

##### Cythonでの高速化1 （Cython1)
コードは変更せずただCythonでコンパイル

```python
% load_ext Cython
%% cython
def find_factor1(n):
    #略
def is_prime1(n):
    #略
```

##### Cythonでの高速化2 （Cython2)
変数を全て型定義

```python
% load_ext Cython
%% cython
def find_factor2(int n):
    cdef int answer = n
    cdef int i
    #略
def is_prime2(n):
    #略
```

##### Cythonでの高速化2-1 （Cython2-1)
for文のiのみを型定義

```python
% load_ext Cython
%% cython
def find_factor2_1(n):
    answer = n
    cdef int i
    #略
def is_prime2_1(n):
    #略
```

##### Cythonでの高速化2-2 （Cython2-2)
answerのみを型定義

```python
% load_ext Cython
%% cython
def find_factor2_2(n):
    cdef int answer = n
    #略
def is_prime2_2(n):
    #略
```

##### Cythonでの高速化2-3 （Cython2-3)
引数nのみを型定義

```python
% load_ext Cython
%% cython
def find_factor2_3(int n):
    #略
def is_prime2_3(int n):
    #略
```

##### Cythonでの高速化2-4 （Cython2-4)
引数n以外を型定義

```python
% load_ext Cython
%% cython
def find_factor2_4(n):
    cdef int answer = n
    cdef int i
    #略
def is_prime2_4(n):
    #略
```

##### Cythonでの高速化3 （Cython3)
findfactor関数をcdef で定義

```python
% load_ext Cython
%% cython
cdef int find_factor3(int n):
    cdef int answer = n
    cdef int i
    #略
def is_prime3(n):
    #略
```

実験は全て素数131071を素数判定するまでの時間を`%timeit`によって計測しました。それぞれ100回実行した処理の平均時間と標準偏差は以下のようになりました。

##### is_prime関数の速度の比較
|        |       Baseline |    Numba    | Cython1<br>コード変更なし    | Cython2   <br> 関数自体も|Cython3  <br> 関数内変数を全て|
|:-----------------:|:------------------:|:------------------:|:------------------:|:------------------:|:------------------:|
| 時間 (ms)            |   8.69±0.2    |        0.49±0.0       |       5.34±0.1        |      0.43±0.0        |      0.43±0.0        |




##### Cythonにて型定義の箇所による速度の比較

|        |     Cython2-1 <br> for文のiのみ   | Cython2-2 <br> answerのみ   | Cython2-3 <br> 引数nのみ |　Cython2-4  <br> 引数n以外|
|:----------------:|:------------------:|:------------------:|:------------------:|:------------------:|
| 時間 (ms)  |   4.81±0.1      |      5.27±0.4       |       6.62±0.1        |      4.85±0.1       |    


Cython2の結果より、Cythonに関しては型情報を全て定義した時一番速くなりました。ここで型定義の箇所の違いによる速度の変化について考察してみます。Cython2-1、Cython2-2の違いからわかるように、変数を型定義した時の効果は前者(イテレータ変数`i`を型定義)の寄与が大きいことが分かります。これはforループで行われているインクリメント計算では演算ごとに型チェックが入る(`int.__add__`を毎回呼び出している処理と同等)ので、多量に実行された場合オーバーヘッドがきいてくるが、型定義によりC言語上で直接`i+1`を計算することでこのオーバーヘッドをなくしているためだと考えられます。

また、Cython2-3に注目すると引数のみの型定義ではCython1と比べて速度が遅くなっています。`find factor`関数を呼ぶ際の引数における動的な型チェックの除去によって速度が速くなるはずと考えていましたが実際には遅くなりました。それらの効果よりも引数の型定義によるコンパイルのオーバーヘッドがあるのかと考え、何も処理しないが引数の型定義の有無のみ違う以下のような関数`hoge1,hoge2`を定義し、速度を比較したところ若干の差が見られました。

```python
%% cython
def hoge1 (n,a,b,c):
    return
def hoge2 (int n, int a, int b, int c):
    return
```

しかし上記の関数でも10nsほどの差であった(1000万ループ回しており標準偏差の値からも差は誤差ではないはず）のに加えて引数が一つの場合はほとんど差がなかったのでこのオーバーヘッドによる差ではないと考えられます。Cython2-4から`find factor`内の変数が全て型定義されている場合は引数を型定義した方が速いことも分かります。以上からCython2-3での遅延の原因は`n % i == 0`の演算において、`n,i`とも型定義されていない場合よりもどちらか片方のみ型定義されている場合の動的チェックによるオーバヘッドの方が大きいからという結論になると思いますが理由はわかりませんでした。

Cython3より`find factor`自体を型定義した結果は、Cython2とそれほど変わりませんでした。`find factor`内の演算処理自体はCython2と同じなので、納得できます。しかしCython3では`find factor`を呼び出すときのオーバーヘッドが除去されるはずなので`find factor`を一度呼ぶだけではほとんど差が現れないが、多量に呼ぶ処理ではCython3の方が速くなるのではないかと考え以下のような関数を作成し検証を行いました。


```python
%% cython

cdef int find_factor3_1(int n):
    #略
def find_factor3_2(int n):
    #略

def is_prime_inf1(n,m):
    for i in range (m):
        find_factor3_1 (n)

def is_prime_inf2(n,m):
    for i in range (m):
        find_factor3_2 (n)
```

##### 一万回find factor を呼び出した時の比較

|        |       find factor3-1(cdef) |    find factor3-2(def)    |
|:-----------------:|:------------------:|:------------------:|
| 時間            |   158 μs |        3.1s     |

想定通り多量に`find factor`を呼ぶ事で有意な差が現れました。よって`find factor`自体を型定義する事で呼び出す際のオーバーヘッドが除去される事が分かりました。



# 実験2 行列積での比較
行列積をベースラインとして速度の高速化を試みました。

#####Pythonによる実装

```python
def dot(a, b):
    n = len(a)
    l = len(a[0])
    m = len(b[0])
    c = np.zeros((n, m))
    for i in range(n):
        for j in range(m):
            for k in range(l):
                c[i][j] += a[i][k] * b[k][j]
    return c
```

#####Numbaによる実装

```python
from numba import jit
@jit
def nm_dot(a, b):
    #略
```

#####Cythonによる実装1 (Cython1)

とりあえず今までの知識を元に型定義を行いcythonで実装。より高速化するために`[i][j]`→`[i,j]`へ変更、型定義はまとめて宣言、`range` → `xrange` へ変更、チェック機能をオフにするなどの工夫を行っています。

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

cython: boundscheck=False
cython: wraparound=False
    
cpdef np.ndarray  c_dot(np.ndarray a, np.ndarray b):
    cdef np.ndarray c
    cdef int n,l,m,i,j,k
    n = len(a)
    l = len(a[0])
    m = len(b[0])
    c = np.zeros((n, m))
    for i in xrange(n):
        for j in xrange(m):
            for k in xrange(l):
                c[i,j] += a[i,k] * b[k,j]
    return c
```

##### Cythonによる実装2 (Cython2)
公式ドキュメントを参考にnumpy配列に対する型定義をより厳密に行いました。

```python
cpdef np.ndarray[dtype=float,ndim = 2]  c_dot2(np.ndarray[dtype=float,ndim = 2] a, np.ndarray[dtype=float,ndim = 2] b):
    cdef np.ndarray[dtype=double,ndim = 2] c
    #略
```

この実験では100*100のランダムNumpy行列を生成し、実験1と同様に`%timeit`で計測しました。処理時間のかかるPythonとCython1の実装では10回、それ以外では100回分の処理の平均と標準偏差は以下のようになりました。



##### 行列積の速度の比較

|        |  Python   |    Python(Numpy)    | Numba    | Cython1 | Cython2 |
|:----------------:|:------------------:|:------------------:|:------------------:|:------------------:|:------------------:|
| 時間 (ms)  |   1220±54      |      0.033       |       1.1       |      729±16  |    2.5      |   

Cython1ではナイーブなPythonでの実装より多少速い程度であるがCython2ではかなり速くなっています。これはCython1での実装では各numpy配列についてnp.ndarrayと型定義した時点ではまだ実際のdtypeや次元が確定せず結局動的チェックの対象となり、Cythonを使わない場合の処理と同等になってしまうためだと考えられます。
結果から分かるように行列積の計算に関してはとにかくNumpyが圧倒的に速かったです。そもそもNumpyを利用したことは多々あっても何故速いのか原理を知らなかったので少し調べたところ(最初にやるべきでした)、Numpyの演算Numpy.adarray型で確定させることで動的チェックを無くし、内部ではC言語による実装がなされているとのことでした。ここでも型付けが出てきて型付けの重要性を痛感しました。
少なくとも行列積などnumpy にモジュールとして実装されている計算においてはそのままnumpyを利用した方がよく、numpy モジュールとして実装されていないかつnumpy配列に実際に何度もアクセスしないと実現できないような場合にCythonを利用するのが良さそうです。





# おわりに
本実験を通して、実用的なCythonの導入では闇雲に型定義するのではなく、Pythonの利点であるコードの柔軟性や可読性、Numpyなどの優良なライブラリの利用を考慮しつつ、効率よくパフォーマンスをあげられる箇所を見定めていくのが良いと感じました。

# 参考資料

- [cython入門](https://qiita.com/en3/items/1f1a609c4d7c8f3066a7)
- [深入りしないCython入門](https://qiita.com/pashango2/items/45cb85390193d97523ca)
- [Cython Document](http://omake.accense.com/static/doc-ja/cython/index.html
8)
- [NumPy Wikipedia](https://ja.wikipedia.org/wiki/NumPy)