2017-07-17: HyからKerasを利用してMNISTを分類する(畳み込みニューラルネット版)

参考

コード

(import [keras [*]]
        [keras.datasets [mnist]]
        [keras.models [Sequential]]
        [keras.layers [Dense Dropout Flatten]]
        [keras.layers [Conv2D MaxPooling2D]]
        [keras [backend :as K]])

(import keras)

(setv batch-size  128
      num-classes 10
      epochs      12)

;; input image dimensions
(setv img-rows 28
      img-cols 28)

;; the data, shuffled and split between train and test sets
(setv mnist-dataset (mnist.load-data)
      x-train (first  (first  mnist-dataset))
      y-train (second (first  mnist-dataset))
      x-test  (first  (second mnist-dataset))
      y-test  (second (second mnist-dataset)))

(if (= (K.image-data-format) "channels_first")
  (setv x-train (.reshape x-train (nth x-train.shape 0) 1 img-rows img-cols)
        x-test  (.reshape x-test  (nth x-test.shape 0)  1 img-rows img-cols)
        input-shape [1 img-rows img-cols])
  (setv x-train (.reshape x-train (nth x-train.shape 0) img-rows img-cols 1)
        x-test  (.reshape x-test  (nth x-test.shape 0)  img-rows img-cols 1)
        input-shape [img-rows img-cols 1]))

;; normalize [0.0 1.0]
(setv x-train (.astype x-train "float32")
      x-test  (.astype x-test  "float32"))
(/= x-train 255)
(/= x-test  255)

(print x-train.shape "train samples")
(print x-test.shape "test samples")

;; one-hot-encoding
(setv y-train (keras.utils.to-categorical y-train num-classes)
      y-test  (keras.utils.to-categorical y-test  num-classes))

;; define model
(defmacro define-sequential [model layers-def compile-options-dict]
  (setv layers
        (list (map (fn (elem) `(.add ~model ~elem)) layers-def)))
  `(do (def ~model (Sequential))
       ~@layers
       (.compile ~model ~@compile-options-dict)
       ~model))

(define-sequential model
  [(Conv2D 32 :kernel-size [3 3] :input-shape input-shape :activation "relu")
   (Conv2D 64 :kernel-size [3 3] :activation "relu")
   (MaxPooling2D :pool-size [2 2])
   (Dropout 0.25)
   (Conv2D 32 :kernel-size [3 3] :activation "relu")
   (Conv2D 64 :kernel-size [3 3] :activation "relu")
   (MaxPooling2D :pool-size [2 2])
   (Dropout 0.25)
   (Flatten)
   (Dense :units 128 :activation "relu")
   (Dropout 0.5)
   (Dense :units num-classes :activation "softmax")]
  {:loss keras.losses.categorical-crossentropy
   :optimizer (keras.optimizers.Adadelta)
   :metrics ["accuracy"]})

(.summary model)

;; run
(def history
  (.fit model x-train y-train
        :batch-size batch-size :epochs epochs :verbose 1 :validation-data [x-test y-test]))

;; test
(def loss-metrics (.evaluate model x-test y-test :batch-size batch-size))
(print "loss-metrics: " loss-metrics)

;; predict
(def classes (.predict-classes model x-test :batch-size batch-size))
(print "predicted classes: " classes)

(def proba (.predict-proba model x-test :batch_size batch-size))
(print "predicted proba: " proba)

50エポックくらい回すと99.55%くらいまでいく。


2017-07-11: cl-random-forestでクロスバリデーション

kaggleに登録してみたのでとりあえずDigit Recognizerをやってみた。このデータセットではテストデータにはラベルが付いてなくて、予測結果を投稿して初めてテストデータに対する正答率が分かる。ので、学習のメタパラメータのチューニングは訓練データでクロスバリデーションすることによって決める。

cl-random-forestにn-foldクロスバリデーションを実装したのでMNISTで試してみる。データセットの分割数を5とすると、

(ql:quickload :cl-random-forest)

;; MNISTのデータを用意する
(defparameter mnist-dim 780)
(defparameter mnist-n-class 10)

(let ((mnist-train (clol.utils:read-data "/home/wiz/tmp/mnist.scale" mnist-dim :multiclass-p t))
      (mnist-test (clol.utils:read-data "/home/wiz/tmp/mnist.scale.t" mnist-dim :multiclass-p t)))

  ;; Add 1 to labels in order to form class-labels beginning from 0
  (dolist (datum mnist-train) (incf (car datum)))
  (dolist (datum mnist-test)  (incf (car datum)))

  (multiple-value-bind (datamat target)
      (clol-dataset->datamatrix/target mnist-train)
    (defparameter mnist-datamatrix datamat)
    (defparameter mnist-target target))
  
  (multiple-value-bind (datamat target)
      (clol-dataset->datamatrix/target mnist-test)
    (defparameter mnist-datamatrix-test datamat)
    (defparameter mnist-target-test target)))
    
;; 並列化を有効化
(setf lparallel:*kernel* (lparallel:make-kernel 4))

;; cross-validation
(defparameter mnist-n-fold 5)

(cross-validation-forest-with-refine-learner
 mnist-n-fold mnist-n-class mnist-dim mnist-datamatrix mnist-target
 :n-tree 500 :bagging-ratio 0.1 :max-depth 10 :n-trial 28 :gamma 10d0 :min-region-samples 5)

これで5-fold cross-validationの平均値が出る。他を固定してn-treebagging-ratioを動かしてみると、 clrf-cv.png 決定木の数n-treeは500くらい以上にしても正答率はほとんど変わらなくなる。論文の通りに、バギング比率bagging-ratioはほとんど正答率に影響を与えないことが分かるので、0.1などの小さい値にすることで計算時間を短縮できる。

ランダムフォレストはあまりメタパラメータに敏感ではないので、割と何も考えずにデフォルト値でやってもうまくいく。逆に言うとあまりチューニングの余地はない。


2017-07-10: JuliaのランダムフォレストライブラリDecisionTree.jlでMNIST

Juliaが速くて機械学習分野で人気と聞いたので試してみた。ランダムフォレストの実装があったのでこれでMNISTデータの分類を試してみる。

インストール

まず必要パッケージをインストールする。

Pkg.add("DecisionTree")
Pkg.add("MNIST")

データの読み込みと前処理

MNISTパッケージにデータが含まれているが、DecisionTreeのデータとして使うにはちょっとした前処理が必要。

# MNISTパッケージの機能でMNISTのデータを読み込む
trainX, trainY = traindata()
testX, testY = testdata()

# 入力行列を転置
trainXt = transpose(trainX)
testXt  = transpose(testX)

# ターゲットが数値だと回帰になってしまうので文字列に直す
function setStringArr(arr,arrstr)
    dataSize = size(arr)[1]
    for i in 1:dataSize
        x=arr[i]
        arrstr[i]="$x"
    end
end

trainYstr = Array{String}(60000)
setStringArr(trainY,trainYstr)

testYstr = Array{String}(10000)
setStringArr(testY,testYstr)

ランダムフォレストのモデル構築

データからランダムフォレストのモデルを作るにはbuild_frest()を使う。引数は、サンプリングする特徴数、木の数、バギング比率、決定木の最大深さ。サンプリングする特徴数は推奨値であるところの元の特徴数784の平方根28を使う。

処理の時間を測るには@timeを頭に付けておけばいいらしい。118秒かかっている。

# using 28 random features, 10 trees, 1.0 portion of samples per tree (optional), and a maximum tree depth of 15 (optional)
@time model = build_forest(trainYstr, trainXt, 28, 10, 1.0, 15)

# 118.467527 seconds (90.35 M allocations: 57.266 GiB, 18.03% gc time)
# Ensemble of Decision Trees
# Trees:      10
# Avg Leaves: 4097.6
# Avg Depth:  15.0

テスト

predictions = apply_forest(model, testXt)
cm = confusion_matrix(testYstr, predictions)

# Classes:  String["0.0", "1.0", "2.0", "3.0", "4.0", "5.0", "6.0", "7.0", "8.0", "9.0"]
# Matrix:   10×10 Array{Int64,2}:
#  963     1    0    0    1    3    8    1    3    0
#    0  1118    3    4    0    1    4    2    3    0
#    7     0  979    5    4    2   11   14    8    2
#    1     0   21  938    1   23    2    8    8    8
#    1     0    5    1  905    1    9    4    8   48
#    5     1    2   37    1  813   13    3    9    8
#    9     3    2    0    4    2  936    0    2    0
#    0     2   19    7    2    0    0  970    5   23
#    5     0   10   18    5   10   11    2  889   24
#    5     4    4   11   11    3    3    7    6  955

# Accuracy: 0.9466
# Kappa:    0.94064143993818

で94.6%くらいの正答率になる。木の数が10にしては時間がかかりすぎている。topを見るかぎり並列化もしていないっぽい。


2017-07-10: HyからKerasを利用してMNISTを分類する

以前HyからChainerを利用してMNISTを分類する例を上げたので、今度はKerasでやってみる。

参考

インストール

pip install tensorflow-gpu
pip install keras

データの用意とモデル定義

(import [keras.models [Sequential]]
        [keras.layers [Dense Dropout Activation]]
        [keras.optimizers [Adam]]
        [keras.datasets [mnist]]
        [keras.utils [np-utils]])

;;; データセット(MNIST)の取得と前処理
(def mnist-dataset (mnist.load-data))
(def X-train (first (first mnist-dataset)))
(def Y-train (second (first mnist-dataset)))
(def X-test (first (second mnist-dataset)))
(def Y-test (second (second mnist-dataset)))

;; 画像を1次元配列化
(setv X-train (.reshape X-train 60000 784)
      X-test  (.reshape X-test  10000 784))

;; 画素を0.0-1.0の範囲に変換
(setv X-train (.astype X-train "float32")
      X-test  (.astype X-test  "float32"))
(/= X-train 255)
(/= X-test  255)

(print X-train.shape "train samples")
(print X-test.shape "test samples")

;; one-hot-encoding
(setv Y-train (np-utils.to-categorical Y-train 10)
      Y-test  (np-utils.to-categorical Y-test  10))

;;; モデル定義
(def model (Sequential))
(.add model (Dense :input-dim 784 :units 512))
(.add model (Activation "relu"))
(.add model (Dropout 0.2))
(.add model (Dense :units 512))
(.add model (Activation "relu"))
(.add model (Dropout 0.2))
(.add model (Dense :units 10))
(.add model (Activation "softmax"))
(.compile model
          :loss "categorical_crossentropy"
          :optimizer (Adam)
          :metrics ["accuracy"])

(.summary model)
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
dense_1 (Dense)              (None, 512)               401920    
_________________________________________________________________
activation_1 (Activation)    (None, 512)               0         
_________________________________________________________________
dropout_1 (Dropout)          (None, 512)               0         
_________________________________________________________________
dense_2 (Dense)              (None, 512)               262656    
_________________________________________________________________
activation_2 (Activation)    (None, 512)               0         
_________________________________________________________________
dropout_2 (Dropout)          (None, 512)               0         
_________________________________________________________________
dense_3 (Dense)              (None, 10)                5130      
_________________________________________________________________
activation_3 (Activation)    (None, 10)                0         
=================================================================
Total params: 669,706
Trainable params: 669,706
Non-trainable params: 0
_________________________________________________________________

マクロを使うとモデル定義を多少すっきりさせられる。

(defmacro define-sequential [model layers-def compile-options-dict]
  (setv layers
        (list (map (fn (elem)
                     `(.add ~model ~elem))
                   layers-def)))
  `(do (def ~model (Sequential))
       ~@layers
       (.compile ~model ~@compile-options-dict)
       ~model))

(define-sequential model
  ;; layer
  [(Dense :input-dim 784 :units 512)
   (Activation "relu")
   (Dropout 0.2)
   (Dense :units 512)
   (Activation "relu")
   (Dropout 0.2)
   (Dense :units 10)
   (Activation "softmax")]
  ;; options
  {:loss "categorical_crossentropy"
   :optimizer (Adam)
   :metrics ["accuracy"]})

訓練、テスト、予測

;; 訓練を実行
(def history (.fit model X-train Y-train :nb-epoch 5 :batch-size 100))

;; テスト
(def loss-metrics (.evaluate model X-test Y-test :batch-size 100))

(for [i (range 4)]
  (.fit model X-train Y-train :nb-epoch 5 :batch-size 100)
  (print "\n" (.evaluate model X-test Y-test :batch-size 100)))

;; 予測(分類結果)
(def classes (.predict-classes model X-test :batch-size 100))
;; 予測(分布)
(def proba (.predict-proba model X-test :batch_size 100))
Epoch 1/10
2017-06-20 13:31:26.460964: W tensorflow/core/platform/cpu_feature_guard.cc:45] The TensorFlow library wasn't compiled to use SSE4.1 instructions, but these are available on your machine and could speed up CPU computations.
2017-06-20 13:31:26.460990: W tensorflow/core/platform/cpu_feature_guard.cc:45] The TensorFlow library wasn't compiled to use SSE4.2 instructions, but these are available on your machine and could speed up CPU computations.
2017-06-20 13:31:26.460997: W tensorflow/core/platform/cpu_feature_guard.cc:45] The TensorFlow library wasn't compiled to use AVX instructions, but these are available on your machine and could speed up CPU computations.
2017-06-20 13:31:26.461001: W tensorflow/core/platform/cpu_feature_guard.cc:45] The TensorFlow library wasn't compiled to use AVX2 instructions, but these are available on your machine and could speed up CPU computations.
2017-06-20 13:31:26.461007: W tensorflow/core/platform/cpu_feature_guard.cc:45] The TensorFlow library wasn't compiled to use FMA instructions, but these are available on your machine and could speed up CPU computations.
2017-06-20 13:31:26.683495: I tensorflow/stream_executor/cuda/cuda_gpu_executor.cc:893] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero
2017-06-20 13:31:26.683701: I tensorflow/core/common_runtime/gpu/gpu_device.cc:940] Found device 0 with properties: 
name: GeForce GTX 750 Ti
major: 5 minor: 0 memoryClockRate (GHz) 1.0845
pciBusID 0000:01:00.0
Total memory: 1.95GiB
Free memory: 1.60GiB
2017-06-20 13:31:26.683719: I tensorflow/core/common_runtime/gpu/gpu_device.cc:961] DMA: 0 
2017-06-20 13:31:26.683724: I tensorflow/core/common_runtime/gpu/gpu_device.cc:971] 0:   Y 
2017-06-20 13:31:26.683733: I tensorflow/core/common_runtime/gpu/gpu_device.cc:1030] Creating TensorFlow device (/gpu:0) -> (device: 0, name: GeForce GTX 750 Ti, pci bus id: 0000:01:00.0)
60000/60000 [==============================] - 76s - loss: 0.2144 - acc: 0.9345    
Epoch 2/10
60000/60000 [==============================] - 6s - loss: 0.1040 - acc: 0.9682     
Epoch 3/10
60000/60000 [==============================] - 6s - loss: 0.0821 - acc: 0.9744     
Epoch 4/10
60000/60000 [==============================] - 6s - loss: 0.0664 - acc: 0.9789     
Epoch 5/10
60000/60000 [==============================] - 6s - loss: 0.0574 - acc: 0.9820     
Epoch 6/10
60000/60000 [==============================] - 6s - loss: 0.0523 - acc: 0.9840     
Epoch 7/10
60000/60000 [==============================] - 6s - loss: 0.0454 - acc: 0.9858     
Epoch 8/10
60000/60000 [==============================] - 6s - loss: 0.0440 - acc: 0.9867     
Epoch 9/10
60000/60000 [==============================] - 6s - loss: 0.0388 - acc: 0.9879     
Epoch 10/10
60000/60000 [==============================] - 6s - loss: 0.0351 - acc: 0.9890     
<keras.callbacks.History object at 0x7f7882ace6d8>
=> (.fit model X-train Y-train :nb-epoch 10 :batch_size 32)
Epoch 1/10
60000/60000 [==============================] - 6s - loss: 0.0341 - acc: 0.9896     
Epoch 2/10
60000/60000 [==============================] - 6s - loss: 0.0350 - acc: 0.9897     
Epoch 3/10
60000/60000 [==============================] - 6s - loss: 0.0332 - acc: 0.9900     
Epoch 4/10
60000/60000 [==============================] - 6s - loss: 0.0326 - acc: 0.9912     
Epoch 5/10
60000/60000 [==============================] - 6s - loss: 0.0325 - acc: 0.9907     
Epoch 6/10
60000/60000 [==============================] - 6s - loss: 0.0293 - acc: 0.9917     
Epoch 7/10
60000/60000 [==============================] - 6s - loss: 0.0277 - acc: 0.9921     
Epoch 8/10
60000/60000 [==============================] - 6s - loss: 0.0294 - acc: 0.9919     
Epoch 9/10
60000/60000 [==============================] - 6s - loss: 0.0257 - acc: 0.9928     
Epoch 10/10
60000/60000 [==============================] - 6s - loss: 0.0291 - acc: 0.9919     
<keras.callbacks.History object at 0x7f7893573a20>

2017-06-06: cl-random-forestでランダムフォレストの決定境界を描いてみる

cl-random-forestでは通常のランダムフォレストに加えて、ランダムフォレストの構造を使って特徴抽出し、それを線形分類器で再学習するという手法を実装している(Global refinement of random forest)。 通常のランダムフォレストに対して、この手法がどういう分類をしているかを見るために、二次元のデータでの実際の分類結果を可視化してみる。

このエントリではXORのデータで、綺麗に分かれている場合かなりオーバーラップしている場合 とでランダムフォレストの決定境界を描いている。データはそれぞれのリンク先にある。

ランダムフォレストを構築

まずはXORのデータからランダムフォレストを構築する。完全なコードはここにある。

(defparameter *target*
  (make-array 100 :element-type 'fixnum
                 :initial-contents
                 '(1 1 0 1 1 0 1 0 1 1 0 0 1 0 0 0 1 0 0 0 0 0 1 1 0 1 1 1 0 0 0 1 1 1 1 1 1 1 0
                   0 1 1 0 0 1 0 0 1 0 1 1 1 1 1 1 0 1 0 1 1 0 0 0 0 1 0 0 1 1 1 1 0 1 0 0 1 0 1
                   0 0 0 0 0 1 0 1 0 1 1 0 0 1 0 0 1 0 0 1 1 0)))

(defparameter *datamatrix*
  (make-array '(100 2) 
              :element-type 'double-float
              :initial-contents
              '((1.7128444910049438d0 -1.1600620746612549d0)
                (1.0938962697982788d0 -1.4356023073196411d0)
                (1.6034027338027954d0 0.6371392607688904d0)
                ...
                (-1.1475433111190796d0 0.9125188589096069d0)
                (0.5465568900108337d0 -0.3704013228416443d0)
                (1.1223702430725098d0 2.0434348583221436d0))))

;;; make random forest
(defparameter *forest*
  (make-forest *n-class* *n-dim* *datamatrix* *target* :n-tree 100 :bagging-ratio 1.0 :max-depth 100))

;;; test random forest
(test-forest *forest* *datamatrix* *target*)
; Accuracy: 100.0%, Correct: 100, Total: 100

可視化

プロットには以前実装したgnuplotのフロントエンドclgplotを使う。(解説記事)

splot-matrix関数でCommon Lispの二次元配列を三次元カラーマップでプロットすることができる。ここではmake-predict-matrixでプロット用の二次元配列を作っている。この関数には行列の一辺のサイズと定義域、さらに学習済みのモデルと訓練データを与える。

まずモデルから定義域全体の予測結果を出して、配列を埋めていく。次に配列中の各訓練データ点に相当する位置にラベルを表す1/-1を入れていく。

;;; plot decision boundary
(ql:quickload :clgplot)

(defun make-predict-matrix (n x0 xn y0 yn forest datamatrix target)
  (let ((x-span (/ (- xn x0) (1- n)))
        (y-span (/ (- yn y0) (1- n)))
        (mesh-predicted (make-array (list n n)))
        (tmp-datamatrix (make-array '(1 2) :element-type 'double-float)))
    ;; mark prediction
    (loop for i from 0 to (1- n)
          for x from x0 by x-span do
            (loop for j from 0 to (1- n)
                  for y from y0 by y-span do
                    (setf (aref tmp-datamatrix 0 0) x
                          (aref tmp-datamatrix 0 1) y)
                    (setf (aref mesh-predicted i j)
                          (if (> (predict-forest forest tmp-datamatrix 0) 0)
                              0.25 -0.25))))
    ;; mark datapoints
    (let* ((range (clgp:seq 0 (1- n)))
           (x-grid (mapcar (lambda (x) (+ (* x x-span) x0)) range))
           (y-grid (mapcar (lambda (y) (+ (* y y-span) y0)) range)))
      (loop for i from 0 to (1- (length target)) do
        (let ((data-i (position-if (lambda (x) (<= (aref datamatrix i 0) x)) x-grid))
              (data-j (position-if (lambda (y) (<= (aref datamatrix i 1) y)) y-grid)))
          (if (and data-i data-j)
              (setf (aref mesh-predicted data-i data-j)
                    (if (> (aref target i) 0) 1 -1))))))
    mesh-predicted))

(defparameter predict-matrix
  (make-predict-matrix 100 -3.5d0 3.5d0 -3.5d0 3.5d0 *forest* *datamatrix* *target*))

(clgp:splot-matrix predict-matrix)

そうするとプロット結果は以下のようになる。

左が綺麗に分かれてるデータで右がオーバーラップしてるデータのときの分類結果を表している。どちらも訓練データを完璧に分類しているが、右の方ではあまりXORっぽくは見えない。特に、左上の外れ値に強く反応してかなりの領域を青側に誤判定してしまっている。

xor-simple.png

次の図は木の数を100から500にし、最大深さを100から15にし、バギング比率を1から0.1にしたときのもの。

xor-shallow-forest.png

上の図の条件下でさらにGlobal refinementを実行したもの。

xor-refine.png

こうして見ると通常のランダムフォレストに対して外れ値に強くなり、より汎化性能が増しているよう見えなくもない。