機械学習アルゴリズムのボトルネックをCythonで改善する話

#この記事について

- Pythonを速くする方法について語ります
- プロファイリングによってCPUバウンドな処理のボトルネックを見つけます
- 見つけたボトルネックをCythonで改善します

# はじめに
先日Bayesian Personalized Ranking (BPR)というレコメンドアルゴリズムを実装しました。
[こちらの論文](http://arxiv.org/abs/1205.2618)の式を参考にコードを書いてみたのですが、実行してみたらちょっと遅すぎて使えなかったため、処理速度の改善に取り組みました。
その時に試したことを備忘録的にまとめます。

# この記事で用いる手法とコード
BPRはユーザ x アイテムの行列の行列分解を与えます。
ユーザ x アイテムの行列$X$を、ユーザ x ファクターの行列$U$とアイテム x ファクターの行列$V$に分解します。

```math
X = U \cdot V^T
```

この問題をどのように解くかはBPRの[元論文](http://arxiv.org/abs/1205.2618)をご覧ください。

この手法を以下のように実装しました。$X$は`scipy.sparse.lil_matrix`で定義しています。分解後の$U$と$V$は`numpy.array`です。
学習に用いるサンプルをブートストラップサンプリングして、その後$U,V$を更新するという流れです。

このコードについて、モデルの学習の部分である`buildModel()`を改善したいと思います。


```py:MFbpr.py
class MFbpr(Recommender):
    '''
    コンストラクタとか他の処理
    '''
    
    def buildModel(self):
        loss_pre = sys.float_info.max
        nonzeros = self.trainMatrix.nnz
        hr_prev = 0.0
        sys.stderr.write("Run for BPR. \n")
        for itr in xrange(self.maxIter):
            start = time.time()
            
            # Each training epoch
            for s in xrange(nonzeros):
                # sample a user
                u = np.random.randint(self.userCount)
                itemList = self.trainMatrix.getrowview(u).rows[0]
                if len(itemList) == 0:
                    continue
                # sample a positive item
                i = random.choice(itemList)
                
                # One SGD step update
                self.update_ui(u, i)
                
    def update_ui(self, u, i):
        # sample a negative item(uniformly random)
        j = np.random.randint(self.itemCount)
        while self.trainMatrix[u, j] != 0:
            j = np.random.randint(self.itemCount)
            
        # BPR update rules
        y_pos = self.predict(u, i)  # target value of positive instance
        y_neg = self.predict(u, j)  # target value of negative instance
        mult = -self.partial_loss(y_pos - y_neg)
        
        for f in xrange(self.factors):
            grad_u = self.V[i, f] - self.V[j, f]
            self.U[u, f] -= self.lr * (mult * grad_u + self.reg * self.U[u, f])
                
            grad = self.U[u, f]
            self.V[i, f] -= self.lr * (mult * grad + self.reg * self.V[i, f])
            self.V[j, f] -= self.lr * (-mult * grad + self.reg * self.V[j, f])
```

```py:exec_bpr.py
from MFbpr import MFbpr

def load_data(ratingFile, testRatio=0.1):
    '''
    データをロードする処理
    '''
    return trainMatrix, testRatings

if __name__ == "__main__":
    # data
    trainMatrix, testRatings = load_data('yelp.rating')

    # settings
    topK = 100
    factors = 64
    maxIter = 10
    maxIterOnline = 1
    lr = 0.01
    adaptive = False
    init_mean = 0.0
    init_stdev = 0.1
    reg = 0.01
    showProgress = False
    showLoss = True

    bpr = MFbpr(trainMatrix, testRatings, topK, factors, maxIter, lr, adaptive, reg, init_mean, init_stdev, showProgress, showLoss)
    bpr.buildModel()

```

なお、コードの全量は[github](https://github.com/tn1031/bpr)で公開しています。

# 実行してみる
とりあえず実行してみます。
###設定
#####データ
- ユーザ数：25,677
- アイテム数：25,815
- レーティング数(評価が存在するサンプルの個数)：627,775

#####ハイパーパラメータ
- ファクター数：64

#####実行環境
VirtualBox上に構築したメモリ2G、プロセッサ2個のubuntuです。

###実行結果
ひとつ目の鍵括弧が1イテレーションにかかった時間[s]、最後の鍵括弧はlossの計算にかかった時間[s]です。

```
> python exex_bpr.py

Data	yelp.rating
#Users	25677, #newUser: 588
#Items	25815
#Ratings	 627775.0 (train), 73167(test), 10056(#newTestRatings)
Run for BPR. 
Iter=0 [128.6936481]	 [-]loss: 624484.357765 [8.16665792465]
Iter=1 [137.202406883]	 [-]loss: 616970.769244 [7.11149096489]
Iter=2 [131.134891987]	 [-]loss: 609361.307524 [7.12593793869]
Iter=3 [134.665620804]	 [-]loss: 601240.628869 [8.45840716362]
Iter=4 [134.722868919]	 [-]loss: 592053.155587 [7.6300880909]
```


約60万サンプルについて計算しているとはいえ、1イテレーション130秒はかかり過ぎな気がします。

# プロファイリング
## 全体のプロファイリング
まずは時間がかかっている部分の特定から始めます。
Pythonプロファイラの`cProfile`を使って処理速度を計測します。
`python -m cProfile -s cumulative ***.py`
と実行すれば時間のかかった処理を降順で表示をしてくれます。


```
> python -m cProfile -s cumulative exex_bpr.py

Data	yelp.rating
#Users	25677, #newUser: 588
#Items	25815
#Ratings	 627775.0 (train), 73167(test), 10056(#newTestRatings)
Run for BPR. 
Iter=0 [244.87996006]	 [-]loss: 624394.802988 [53.4806399345]
Iter=1 [248.624686956]	 [-]loss: 616876.50976 [48.6073460579]
Iter=2 [253.822627068]	 [-]loss: 609269.820414 [52.5446169376]
Iter=3 [261.039247036]	 [-]loss: 601207.904104 [53.8221797943]
Iter=4 [260.285779953]	 [-]loss: 592212.148141 [54.9556028843]
         369374621 function calls (368041918 primitive calls) in 1808.492 seconds

   Ordered by: cumulative time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.171    0.171 1808.493 1808.493 exex_bpr.py:3(<module>)
        1   34.307   34.307 1532.124 1532.124 MFbpr.py:40(buildModel)
  3067088  535.383    0.000  830.046    0.000 MFbpr.py:92(update_ui)
  6209078   48.829    0.000  433.564    0.000 lil.py:228(__getitem__)
  3292937   26.140    0.000  376.631    0.000 lil.py:191(getrowview)
  3292939   66.488    0.000  346.337    0.000 lil.py:85(__init__)
        5    0.000    0.000  263.411   52.682 MFbpr.py:117(_showLoss)
        5   22.098    4.420  263.410   52.682 MFbpr.py:128(loss)
        1    9.443    9.443  242.550  242.550 exex_bpr.py:9(load_data)
```

`Ordered by: cumulative time`という行の下から、関数名と処理に要した時間が羅列されています。これを見ると`update_ui`という関数が3,067,088回呼ばれてトータル535.383秒かかっていることがわかります。
(まあ、`buildModel`の中ではほとんど`update_ui`しか呼んでないので当然といえば当然ですが．．．)

なお、1イテレーションの実行時間が先程より大きくなっているのはcProfileのオーバーヘッドです。


## 行単位のプロファイリング
[line_profiler](https://pypi.python.org/pypi/line_profiler/)を使えば注目している関数について行単位で計測できます。
pipでインストールできます。

```
> pip install line_profiler
```

`line_profiler`でプロファイルするためには注目している関数に`@profile`のデコレータを付ける必要があります。

```py:MFbpr.py
    @profile
    def update_ui(self, u, i):
        # sample a negative item(uniformly random)
        j = np.random.randint(self.itemCount)
        while self.trainMatrix[u, j] != 0:
            j = np.random.randint(self.itemCount)

        # BPR update rules
        y_pos = self.predict(u, i)  # target value of positive instance
        y_neg = self.predict(u, j)  # target value of negative instance
        mult = -self.partial_loss(y_pos - y_neg)

        for f in xrange(self.factors):
            grad_u = self.V[i, f] - self.V[j, f]
            self.U[u, f] -= self.lr * (mult * grad_u + self.reg * self.U[u, f])

            grad = self.U[u, f]
            self.V[i, f] -= self.lr * (mult * grad + self.reg * self.V[i, f])
            self.V[j, f] -= self.lr * (-mult * grad + self.reg * self.V[j, f])

```

あとは`kernprof`コマンドで実行するだけです。

```
> kernprof -l -v exex_bpr.py

Data	yelp.rating
#Users	25677, #newUser: 588
#Items	25815
#Ratings	 627775.0 (train), 73167(test), 10056(#newTestRatings)
Run for BPR. 
Iter=0 [953.993126154]	 [-]loss: 624406.940531 [7.50253987312]
Iter=1 [962.82383585]	 [-]loss: 616855.373858 [7.96375918388]
Wrote profile results to exex_bpr.py.lprof
Timer unit: 1e-06 s

Total time: 1082.55 s
File: MFbpr.py
Function: update_ui at line 92

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    92                                               @profile
    93                                               def update_ui(self, u, i):
    94                                                   # sample a negative item(uniformly random)
    95   1226961      8241361      6.7      0.8          j = np.random.randint(self.itemCount)
    96   1228124     39557350     32.2      3.7          while self.trainMatrix[u, j] != 0:
    97      1163         6123      5.3      0.0              j = np.random.randint(self.itemCount)
    98                                                       
    99                                                   # BPR update rules
   100   1226961      9495381      7.7      0.9          y_pos = self.predict(u, i)  # target value of positive instance
   101   1226961      4331378      3.5      0.4          y_neg = self.predict(u, j)  # target value of negative instance
   102   1226961     10002355      8.2      0.9          mult = -self.partial_loss(y_pos - y_neg)
   103                                                   
   104  79752465    126523856      1.6     11.7          for f in xrange(self.factors):
   105  78525504    161882071      2.1     15.0              grad_u = self.V[i, f] - self.V[j, f]
   106  78525504    191293505      2.4     17.7              self.U[u, f] -= self.lr * (mult * grad_u + self.reg * self.U[u, f])
   107                                                           
   108  78525504    137871145      1.8     12.7              grad = self.U[u, f]
   109  78525504    194033291      2.5     17.9              self.V[i, f] -= self.lr * (mult * grad + self.reg * self.V[i, f])
   110  78525504    199315454      2.5     18.4              self.V[j, f] -= self.lr * (-mult * grad + self.reg * self.V[j, f])

```

`% Time`のカラムを見ると、`update_ui`の処理時間のうち、**90%以上がfor文以下で費やされている**ことがわかりました。

1イテレーションが1,000秒弱と、オーバーヘッドが大変なことになっています。

#Cythonで高速化
##実装
先ほどのforループをCythonで書き直します。
Pythonが遅い原因のひとつが、Pythonは動的型付けであるということです。変数を参照するたびに型をチェックしているので、変数の参照が多いプログラムはこの操作にかかる時間が無視できなくなってきます。

Cythonを使って書けばc言語のように変数の型を定義できます。いちいち型の確認が行われないのでかなりの高速化が期待できます。

変数の宣言は`cdef`で行います。計算に用いるすべての変数の型を定義しておきます。

```py:fastupdfn.pyx
import numpy as np
cimport numpy as np
cimport cython

DOUBLE = np.float64
ctypedef np.float64_t DOUBLE_t

cpdef c_upd(np.ndarray[DOUBLE_t, ndim=2] U, 
          np.ndarray[DOUBLE_t, ndim=2] V,
          double mult,
          double lr,
          double reg,
          unsigned int u,
          unsigned int i,
          unsigned int j,
          unsigned int factors):
    cdef unsigned int f
    cdef double grad_u, grad
    for f in xrange(factors):
        grad_u = V[i, f] - V[j, f]
        U[u, f] -= lr * (mult * grad_u + reg * U[u, f])
        
        grad = U[u, f]
        V[i, f] -= lr * (mult * grad + reg * V[i, f])
        V[j, f] -= lr * (-mult * grad + reg * V[j, f])
        
    return U, V
```

呼び出し元は以下のように書き換えます。

```py:MFbpr.py
import fastupd    # 追加

class MFbpr(Recommender):
    """
    中略
    """
    def update_ui(self, u, i):
        # sample a negative item(uniformly random)
        j = np.random.randint(self.itemCount)
        while self.trainMatrix[u, j] != 0:
            j = np.random.randint(self.itemCount)
            
        # BPR update rules
        y_pos = self.predict(u, i)  # target value of positive instance
        y_neg = self.predict(u, j)  # target value of negative instance
        mult = -self.partial_loss(y_pos - y_neg)
       
        # Cython実装を呼び出す
        self.U, self.V = fastupd.c_upd(self.U, self.V, mult, self.lr, self.reg, u, i, j, self.factors) 

```

また、Cython実装をコンパイルするための`setup.py`も用意する必要があります。

```py:setup.py
from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext

setup(
    cmdclass={'build_ext': build_ext},
    ext_modules=[Extension("fastupd", ["fastupdfn.pyx"])]
)
```

準備が済んだらコンパイルです。次のコマンドでコンパイルします。

```
> python setup.py build_ext --inplace
```

##実行
実行します。

```
> python exec_bpr.py

Data	yelp.rating
#Users	25677, #newUser: 588
#Items	25815
#Ratings	 627775.0 (train), 73167(test), 10056(#newTestRatings)
Run for BPR. 
Iter=0 [38.7308020592]	 [-]loss: 624282.459815 [8.64863014221]
Iter=1 [36.7822458744]	 [-]loss: 616740.942017 [8.62252616882]
Iter=2 [35.8996829987]	 [-]loss: 609119.520253 [7.9108710289]
Iter=3 [35.1236720085]	 [-]loss: 600992.740207 [8.14179396629]
Iter=4 [34.9632918835]	 [-]loss: 591877.909123 [8.81325411797]
```

1イテレーションの実行時間が40秒かからない程度になりました。最初と比べると**3~4倍速くなった**ことになります。

`line_profiler`の結果も載せておきます。

```
> kernprof -l -v exex_bpr.py

Data	yelp.rating
#Users	25677, #newUser: 588
#Items	25815
#Ratings	 627775.0 (train), 73167(test), 10056(#newTestRatings)
Run for BPR. 
Iter=0 [66.7851500511]	 [-]loss: 624400.680806 [7.92221903801]
Iter=1 [62.5339269638]	 [-]loss: 616866.244211 [7.85720801353]
Iter=2 [65.6408250332]	 [-]loss: 609235.226048 [8.48338794708]
Iter=3 [66.0613160133]	 [-]loss: 601140.035318 [8.52119803429]
Iter=4 [66.5882000923]	 [-]loss: 592026.927719 [8.32318806648]
Wrote profile results to exex_bpr.py.lprof
Timer unit: 1e-06 s

Total time: 164.139 s
File: MFbpr.py
Function: update_ui at line 93

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    93                                               @profile
    94                                               def update_ui(self, u, i):
    95                                                   # sample a negative item(uniformly random)
    96   3066856     17642519      5.8     10.7          j = np.random.randint(self.itemCount)
    97   3069840     79530375     25.9     48.5          while self.trainMatrix[u, j] != 0:
    98      2984        15027      5.0      0.0              j = np.random.randint(self.itemCount)
    99                                                       
   100                                                   # BPR update rules
   101   3066856     17651846      5.8     10.8          y_pos = self.predict(u, i)  # target value of positive instance
   102   3066856     10985781      3.6      6.7          y_neg = self.predict(u, j)  # target value of negative instance
   103   3066856     18993291      6.2     11.6          mult = -self.partial_loss(y_pos - y_neg)
   104                                                  
   105   3066856     19320147      6.3     11.8          self.U, self.V = fastupd.c_upd(self.U, self.V, mult, self.lr, self.reg, u, i, j, self.factors) 

```

もともと90%以上を費やしていた行列更新の部分が11.8%に減少しました。Cythonの効果が確認できます。

#まとめ
`cPrifile`や`line_profiler`を用いたプロファイリングによって処理のボトルネックが発見でき、`Cython`を使って改善しました。
1箇所書き換えただけですが、4倍弱速くなりました。

ちなみに、Cythonを適用したコードは[github](https://github.com/tn1031/bpr)のw_cythonブランチに置きました。そのうちmasterにマージするかもしれません。

#おまけ？
for文の部分は、行列の要素を独立に更新しているので完全に並列実行可能です。Cythonには`prange()`という、for文の処理をマルチスレッドで実行してくれる関数があるので、これを適用することでもう少し改善できるかもしれません。