Common Lispによる深層学習ライブラリMGL

昨今の機械学習界隈では何故かPythonやRのライブラリが多いのだが、Common Lispも高速なネイティブコンパイラを持ち、C/C++のライブラリも呼べるので、機械学習にも問題なく使える。
実際にCommon Lispで実装された深層学習を可能とするライブラリにMGLがあり、GithubでMITライセンスで公開されている。そんなわけで、先月のLispmeetup #39でMGLを紹介してきた。これがそのときのスライド。

Lispmeetup #39 MGLの紹介： Common Lispによるディープラーニング from Satoshi imai

MGLの作者Gábor Melis氏はKaggleの2014年の機械学習コンペHiggs Boson Machine Learning Challengeの優勝者で、今はDeepMindにいるらしい。このコンペのコードが公開されているが、MGLが使われていることが分かる。またMGLには深層学習以外にも進化計算や自然言語処理などのパッケージもあるらしい。ドキュメントもちゃんと書かれており、利用するのに必要な基本的材料は揃っているように見える。
MGLは行列演算に以前紹介したMGL-MATを使う。よってCPUであればOpenBLASã‚„MKLなどを、GPUであればCUDAを利用して高速に学習できる。以前、自分でも多層ニューラルネットを実装してみたが、ごくごくナイーブな実装なので実用とは言えなかった。それに対して、MGLはPython/C++のライブラリと比べても遜色ない性能を示す（後述のベンチマーク参照）。
なお、MGLは少なくともLinuxとMacで動くはずだが、自分はLinuxでしか試してないので、Macでも動いたという人は教えてほしい。*1
GPUは別に必須というわけではなく、CPUのみで使う場合はCUDAのインストールも不要である。

• 現状MGLでできること
  • 多層ニューラルネットワーク(ReLU、Maxoutなど各種活性化関数、Dropout完備)
  • リカレントニューラルネットワーク、LSTM
  • 制約ボルツマンマシン(RBM)ベースの教師なし事前学習(DBN、DBM)
  • cl-cudaによるGPU計算、とてもはやい
  • LLA経由でOpenBLAS等を使ってCPU計算、かなりはやい
• できないこと
  • マルチGPUはできないっぽい
  • 畳み込みニューラルネットワーク*2

MGLのインストール

MGL-MATはLLAとcl-cudaに依存する。LLAとMGL-MATのインストールについては以前の記事を参照。

• Common Lispで行列演算: LLA(Lisp Linear Algebra)を使う
• Common Lispで行列演算(2): MGL-MATを使う

MGL-MATのインストールさえできていればMGLのインストールは簡単で、GithubからMGLのソースをダウンロードしてきてquickloadするだけである。その他の依存ライブラリはQuicklispから入る。

cd ~/quicklisp/local-projects
git clone https://github.com/melisgl/mgl.git

(ql:quickload :mgl)

MGLに付属のMNISTのサンプルが長いので整理したものを作った。以下ではこのサンプルコードに従って話を進めていく。

cd ~/quicklisp/local-projects
git clone https://github.com/masatoi/mgl-user.git

(ql:quickload :mgl-user)

この中の src/example.lisp を読むと大体の流れが分かるかと思う。

MGL-MATの設定

深層学習では数値の精度はそれほど必要ではなく、単精度で十分らしい。*3 そこでMGL-MATで扱う数値のデフォルトの型を指定する変数*default-mat-ctype*の値を:doubleから:floatにしておく。特にGPU計算ではこの設定の有無で速度に大きな違いが出る。また、CUDAを使うかどうかを変数*cuda-enabled*によって制御できる。

;;; GPU向けの設定
(setf *default-mat-ctype* :float)
(setf *cuda-enabled* t)

MNISTデータの読み込み

MNISTデータの配布ページから4つのファイルをダウンロードしてきて、適当なディレクトリにgunzipなどで展開する。
次にそのディレクトリへのパスを設定する。

(defparameter *mnist-data-dir* "/path/to/mnist-data/")

この状態で以下を実行すると、変数*training-data*と*test-data*にデータが格納される。

(progn (training-data) (test-data) 'done)

これらのデータ読み込み関数はload-mnist.lispで定義されている。当然ながら他のフォーマットのデータを扱うときはこの部分は自分で書く必要がある。
さて、読み込んだデータはデータ点を表す構造体datumのベクタになっている。

MGL-USER> (aref *training-data* 0)
#S(DATUM
   :ID 1
   :LABEL 5
   :ARRAY #<MAT
            784 B #(0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
                    ...
                    0.0 0.0 0.011764706 0.07058824 0.07058824 0.07058824
                    0.49411765 0.53333336 0.6862745 0.101960786 0.6509804 1.0
                    0.96862745 0.49803922 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
                    0.0 0.0 0.0 0.11764706 0.14117648 0.36862746 0.6039216
                    ...
                    0.0 0.0 0.0)>)

datumのスロットは通し番号idと、正解ラベルlabel、そして実際のデータarrayの3つ。arrayが普通のCommon Lispの配列ではなくて、MGL-MATのMAT構造体になっていることに注意する。

モデル定義

フィードフォワードニューラルネットワークはfnnクラスのオブジェクトとして作られる。fnnのmake-instanceをラップするマクロとしてbuild-fnnがあるので、それを使って隠れ層2層のネットワークを作ってみるとこうなる。

(defparameter my-fnn
  (build-fnn (:class 'fnn :max-n-stripes 100)                     ; ミニバッチのサイズ100
    (inputs     (->input :size 784))                              ; 入力次元数784
    (f1         (->relu (->activation inputs :size 256)))         ; ユニット数256、活性化関数はReLU
    (f2         (->relu (->activation f1     :size 256)))         ; ユニット数256、活性化関数はReLU
    (prediction (->softmax-xe-loss (->activation f2 :size 10))))) ; 出力は10通りの値、活性化関数はソフトマックス関数

ここで定義されたネットワークの中身を見てみると、以下のようになっている。

MGL-USER> (describe my-fnn)
#<FNN {1014AC5903}>
BPN description:
  CLUMPS = #(#<->INPUT INPUTS :SIZE 784 1/100 :NORM 0.00000>
             #<->ACTIVATION (#:G1298 :ACTIVATION) :STRIPES 1/100 :CLUMPS 4>
             #<->RELU F1 :SIZE 256 1/100 :NORM 0.00000>
             #<->ACTIVATION (#:G1299 :ACTIVATION) :STRIPES 1/100 :CLUMPS 4>
             #<->RELU F2 :SIZE 256 1/100 :NORM 0.00000>
             #<->ACTIVATION (#:G1300 :ACTIVATION) :STRIPES 1/100 :CLUMPS 4>
             #<->SOFTMAX-XE-LOSS PREDICTION :SIZE 10 1/100 :NORM 0.00000>)
  N-STRIPES = 1
  MAX-N-STRIPES = 100

ネットワークの実体はlumpという、レイヤーのような役割を持つクラスのオブジェクトの集合で、lumpには->が接頭辞として付いている。またlumpは同名のlumpを生成する関数を持っている。例えば (->input :size 784) はサイズ784の入力層を表す->inputクラスのインスタンスを作る関数を呼ぶ。
build-fnnのキーワードオプション:classは生成するネットワークのオブジェクトのクラスで、:max-n-stripesがミニバッチのバッチサイズを表す。build-fnnの本体部分はlet*みたいなもので(実際let*に展開される)、lumpにつける名前とそのlumpを生成する関数呼び出しの対が並んでいる構造になっており、let*と同様に以前に束縛された変数を後の定義中で使うことができる。
lumpには->inputや->reluなどのユニットの値を持つものと、層間の線型写像を担当する->activationがある。この例では->input、->relu、->softmax-xe-lossがユニットの値を持っている。具体的には、順伝搬の値が入るnodesというスロットと、逆伝搬の値が入るderivativesというスロットを持っており、どちらもバッチサイズ×ユニット数の行列になっている。
例えば入力層のlumpを表示してみると以下のようになっており、nodesとderivativesに784次元の入力データをバッチサイズ分束ねて行列としたものが入っていることが分かる。

MGL-USER> (describe (aref (clumps my-fnn) 0))

#<->INPUT INPUTS :SIZE 784 100/100 :NORM 91.96476>
  [standard-object]

Slots with :INSTANCE allocation:
  NAME               = INPUTS
  SIZE               = 784
  NODES              = #<MAT 100x784 AF #2A((0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0..
  DERIVATIVES        = #<MAT 100x784 A #2A((-1.6850416e-7 1.7476857e-6 3.683108e-6..
  DEFAULT-VALUE      = 0
  SHARED-WITH-CLUMP  = NIL
  X                  = #<->INPUT INPUTS :SIZE 784 100/100 :NORM 91.96476>
  DROPOUT            = NIL
  MASK               = NIL

さて、次は隠れ層についてだが、隠れ層のlumpはその層における活性化関数の種類に応じて->relu、->sigmoid、->tanh、->maxoutなどの中から選ぶ。最近は深層学習にはReLUやMaxoutを使うことが多い。
出力層のlumpは、MNISTが分類問題であるので、活性化関数をソフトマックス関数、損失関数をクロスエントロピーとする ->softmax-xe-loss を使っている。->softmax-xe-lossにはtargetというスロットがあり、正解ラベルをバッチサイズ分並べたリストが入る。これが回帰問題であれば、活性化関数は恒等写像で損失関数は二乗誤差となり、代わりに ->squared-difference を使うことになる。回帰については別の記事で紹介できればと思う。

• >activationは前の層の出力から次の層の入力であるユニット活性を計算するためのlumpで、これ自体がlumpの集合を持つサブネットワークになっている。ここに学習対象である層間の重みやバイアスが入っている。

MaxoutとDropout

活性化関数にMaxoutを使うとネットワークを多層にしたときでも学習がうまくいくことが知られている。Maxoutは複数の線形ユニットからなる関数で、複数の直線からその入力における最大値を選ぶことによって下に凸な関数になる。そのため何個の直線を使うかをパラメータ:group-sizeで指定してやる必要がある。例えば5本の直線を使うなら以下のように書ける。

(->max activations :group-size 5)

Dropoutは学習のたびにユニットの一部を一定確率で無効化する。これは沢山のサブネットワークでアンサンブル学習しているようなもので、強い正則化の効果がある。例えば上の活性化関数でユニットの半分をランダムに無効化するときは、このように書ける。

(->dropout (->max activations :group-size 5) :dropout 0.5)

これらをまとめると、ネットワーク全体はこのように書ける。

(defparameter fnn-maxout-dropout
  (let ((group-size 5))
    (build-fnn (:class 'fnn :max-n-stripes 100)
      (inputs (->input :size 784 :dropout 0.2))                    ; 入力層は0.2の確率で無効化
      (f1-activations (->activation inputs :name 'f1 :size 1200))
      (f1* (->max f1-activations :group-size group-size))
      (f1 (->dropout f1* :dropout 0.5))                            ; 隠れ層は0.5の確率で無効化
      (f2-activations (->activation f1 :name 'f2 :size 1200))
      (f2* (->max f2-activations :group-size group-size))
      (f2 (->dropout f2* :dropout 0.5))
      (f3-activations (->activation f2 :name 'f3 :size 1200))
      (f3* (->max f3-activations :group-size group-size))
      (f3 (->dropout f3* :dropout 0.5))
      (prediction (->softmax-xe-loss (->activation f3 :name 'prediction :size 10)
                                     :name 'prediction)))))

ネットワークへの入出力の設定

こうしてネットワークのオブジェクトを生成したわけだが、学習や予測を行なうためには、まずこのネットワークに入出力を設定してやる必要がある。
入力は->inputのnodesスロットに、出力は->softmax-xe-lossのtargetスロットにそれぞれ設定する。これを行うメソッドがset-inputであり、例えば以下のように定義される。

(defmethod set-input (samples (bpn fnn))
  (let* ((inputs     (find-clump 'inputs     bpn))
         (prediction (find-clump 'prediction bpn)))
    (clamp-data samples (nodes inputs))                       ; バッチサイズ分の入力を設定
    (setf (target prediction) (label-target-list samples))))  ; バッチサイズ分の出力(正解ラベル)を設定

set-input によってネットワークに正しくデータの値が設定されるのであれば、データの表現は別に何でもいいということが分かる。

学習部分の定義

次に、学習を行なう関数を定義する。
ここでどんなアルゴリズムを使うか、学習の途中経過をモニタリングするか等を決める。例えばモーメンタムSGDを使って最適化するときは、以下のように書ける。

(defun train-fnn (fnn training &key
                                 (n-epochs 3000)
                                 (learning-rate 0.1) (momentum 0.9))
  (let ((optimizer (make-instance 'segmented-gd-optimizer
                      :segmenter
                      (constantly
                       (make-instance 'sgd-optimizer
                          :learning-rate learning-rate
                          :momentum momentum
                          :batch-size (max-n-stripes fnn)))))
        (learner (make-instance 'bp-learner :bpn fnn))
        (dateset (make-sampler training :n-epochs n-epochs)))
    (minimize optimizer learner :dataset dateset)
    fnn))

最適化の本体はminimizeで、最適化アルゴリズムを表すoptimizerと学習主体を表すlearnerを引数に取る。ここではoptimizerのところはモーメンタムSGDを指定しているが、他にもADAMなども使える。最適化アルゴリズムは学習率やモーメンタムなど独自のパラメータを持つ。
optimizerやlearnerには学習の途中経過を出力するモニター関数を噛ませることができて、例えば訓練データとテストデータに対する正答率を表示する関数log-bpn-test-errorを最適化の途中で一定間隔で実行する場合にはこのようにする。

(defun train-fnn-with-monitor (fnn training test
                               &key
                                 (n-epochs 3000)
                                 (learning-rate 0.1) (momentum 0.9))
  (let ((optimizer (monitor-optimization-periodically
                    (make-instance 'segmented-gd-optimizer-with-data
                       :training training
                       :test test
                       :segmenter (constantly
                                   (make-instance 'sgd-optimizer
                                      :learning-rate learning-rate
                                      :momentum momentum
                                      :batch-size (max-n-stripes fnn))))
                    '((:fn log-bpn-test-error :period log-test-period)
                      (:fn reset-optimization-monitors
                       :period log-training-period
                       :last-eval 0))))
        (learner (make-instance 'bp-learner :bpn fnn))
        (dateset (make-sampler training :n-epochs n-epochs)))
    (minimize optimizer learner :dataset dateset)
    fnn))

これらを使って、学習の進行過程全体はこのように書ける。

(defun train-fnn-process (fnn training test &key (n-epochs 30) (learning-rate 0.1) (momentum 0.9))
  (with-cuda* ()
    (repeatably ()
      (init-bpn-weights fnn :stddev 0.01) ; fnnの重みを初期化する
      (train-fnn-with-monitor
       fnn training test :n-epochs n-epochs :learning-rate learning-rate :momentum momentum)))
  fnn)

ここでwith-cuda*はMGL-MATの記事でも出てきたCUDAによる計算を有効化するためのマクロで、init-bpn-weightsは正規分布に従って重みの初期値を設定する関数。n-ephochは訓練データを何周するかを表す。

学習プロセスを実行

先に定義したmy-fnnを実際に学習させるとmy-fnn中の重みパラメータが破壊的に更新される。学習プロセス全体の出力は、

MGL-USER> (train-fnn-process my-fnn *training-data* *test-data* :n-epochs 10)
...
2016-05-22 18:40:16: ---------------------------------------------------
2016-05-22 18:40:17: training at n-instances: 600000
2016-05-22 18:40:17: test at n-instances: 600000
2016-05-22 18:40:18: pred. train bpn PREDICTION acc.: 99.75% (10000)
2016-05-22 18:40:18: pred. train bpn PREDICTION xent: 9.027d-5 (10000)
2016-05-22 18:40:18: pred. test  bpn PREDICTION acc.: 98.10% (10000)
2016-05-22 18:40:18: pred. test  bpn PREDICTION xent: 8.772d-4 (10000)
...

となって最終的にテストデータに対して98.10%の正答率であることが分かる。
Dropoutを入れた以下のネットワークで100エポックほど回すと98.96%までいった。

(defparameter fnn-relu-dropout
  (build-fnn (:class 'fnn :max-n-stripes 100)
    (inputs (->input :size 784 :dropout 0.2))
    (f1-activations (->activation inputs :name 'f1 :size 1200))
    (f1* (->relu f1-activations))
    (f1 (->dropout f1* :dropout 0.5))
    (f2-activations (->activation f1 :name 'f2 :size 1200))
    (f2* (->relu f2-activations))
    (f2 (->dropout f2* :dropout 0.5))
    (f3-activations (->activation f2 :name 'f3 :size 1200))
    (f3* (->relu f3-activations))
    (f3 (->dropout f3* :dropout 0.5))
    (prediction (->softmax-xe-loss (->activation f3 :name 'prediction :size 10)
                                   :name 'prediction))))

サンプルファイルのコメントによると、DBMで事前学習してからDropoutありでファインチューニングすることによって99.22%までいくらしい。

予測

データセットの検証にはmonitor-bpn-resultsという関数が用意されており、例えば以下のように使う。

(defun test-fnn (fnn test)
  (monitor-bpn-results (make-sampler test :max-n (length test))
                       fnn
                       (make-datum-label-monitors fnn)))

MGL-USER> (test-fnn my-fnn *test-data*)
(#<CLASSIFICATION-ACCURACY-COUNTER bpn PREDICTION acc.: 98.10% (10000)>
 #<CROSS-ENTROPY-COUNTER bpn PREDICTION xent: 8.772d-4 (10000)>)

個々のデータを予測するには、ネットワークの入力層のnodesに値を設定し、forward関数をネットワークを引数として呼んだ後、出力層のnodesを見れば予測結果が入っている。

(defun predict-datum (fnn datum)
  (let* ((a (datum-array datum))
         (len (mat-dimension a 0))
         (input-nodes (nodes (find-clump 'inputs fnn)))
         (output-nodes (nodes (find-clump 'prediction fnn))))
    ;; set input
    (loop for i from 0 to (1- len) do
      (setf (mref input-nodes 0 i) (mref a i)))
    ;; run
    (forward fnn)
    ;; return output
    (reshape output-nodes (mat-dimension output-nodes 1))))

例えば、訓練データの最初のデータを学習済みのmy-fnnによって予測してみると、以下のような合計が1になるような確率リストが表示される。

MGL-USER> (predict-datum my-fnn (aref *training-data* 0))
#<MAT 0+10+990 - #(1.6357367e-20 1.007323e-17 1.829405e-18 1.9152902e-5
                   3.730603e-23 0.9999808 4.772671e-19 6.1250017e-20
                   1.6942288e-18 2.232391e-16)>

このうち最も確率が高いクラス5が予測結果となり、データの正解ラベルも5なので正解していることが分かる。

ベンチマーク

MNISTを多層ニューラルネットで学習する例をTensorflow、Chainerと比較してみる。
LLAでOpenBLASを使っているので、NumpyでもOpenBLASを使うようにする。(参考: http://qiita.com/unnonouno/items/8ab453a1868d77a93679)
~/.numpy-site.cfg に

[openblas]
libraries = openblas
library_dirs = /usr/lib/openblas-base/
include_dirs = /usr/include/
runtime_library_dirs = /usr/lib/openblas-base/

と書いて pip install numpy するとOpenBLASを使うようになる。
Tensorflowの多層ニューラルネットの例はこれを元にした。ChainerにはMNISTのサンプルが付属するのでそれを元にした。
隠れ層は2層。どちらもユニット数は256、活性化関数をReLU、最適化アルゴリズムをモーメンタムSGDとし、バッチサイズ100で15エポック学習したときの、学習部分のみに要した時間で比較する。環境は Core i5 4670、GeForce GTX 750Ti。
結果はこのようになる。

まとめ

• MGLの特徴
  • 速い
  • 行列演算に関わる部分以外はすべてCommon Lispで実装されており、Common Lispで拡張可能
  • 実行速度を出すために副作用を使いまくっているので関数型的とはいえない
• 個々のケースにおいて実装しなければならないこと。決めなければならないこと
  • データを読み込んで何らかのデータ構造に入れる部分
  • バッチサイズ分のデータをネットワークに設定する部分 (set-inputメソッド)
  • ネットワークの構造
    • 隠れ層の層数、各層のユニット数、活性化関数の選択、Maxoutの場合はgroup-size
    • Dropoutをやる場合はユニットを無効化する確率
  • Optimizerにおけるアルゴリズムの選択とパラメータの指定
    • 例えばmomentum SGDなら学習率とモーメンタムを決める
  • モニタリング関数、ロギング関数など

TODO

• MNISTのサンプルではさらにL2正則化をやっているのでそこを詳しく見ていく。
• 回帰問題の例
• RNN/LSTMの例
• Dropconnectの実装
• 前処理で精度が上がるか: Grobal contrast normalizationなど http://www.slideshare.net/Takayosi/miru2014-tutorial-deeplearning-37219713 P.53

*1:ちなみに当方の環境はハードがCore i5 4670、GeForce GTX 750Tiの自作機、ソフトがUbuntu 14.04(64bit)、CUDA 7.5、OpenBLAS 0.2.8、SBCL 1.3.1となっている。

*2:とはいえ、最大プーリングはMaxoutとほとんど同じ効果という話もあるので、Maxout+DropoutのDNNでも平行移動に対する耐性はある

*3:最近のGPUでは半精度でSIMD演算することでさらに高速化しているらしい