2017-05-17: ChainerのMNISTの例をHyに翻訳してみる

前回Hyのチュートリアルをやってみたので、今度はChainerのサンプルコードをHyに翻訳してみる。

ほとんど逐語訳で済むが、注意点としては、

  • アンダーバーはハイフンに直す
  • オブジェクトの属性はそのままobj.attrでアクセスできる
  • メソッド呼び出しは(obj.method arg1 arg2)でも(.method obj arg1 arg2)でもよい
  • Hyでの多重代入(分配束縛)のやりかたはよく分からなかった。マクロを書けば作れる

コード

;; -*- coding:utf-8; mode:hy -*-

(try
 (import matplotlib)
 (except [ImportError]))

(import [argparse]
        [chainer]
        [chainer.functions :as F]
        [chainer.links :as L]
        [chainer [training]]
        [chainer.training [extensions]])

(defclass MLP [chainer.Chain]
  (defn --init-- [self n-units n-out]
    (.--init-- (super MLP self)
               :l1 (L.Linear None n-units)
               :l2 (L.Linear None n-units)
               :l3 (L.Linear None n-out)))
  (defn --call-- [self x]
    (setv h1 (F.relu (self.l1 x))
          h2 (F.relu (self.l2 h1)))
    (self.l3 h2)))

(defn main []
  (setv parser (argparse.ArgumentParser :description "Chainer example: MNIST"))
  (.add-argument parser "--batchsize" "-b"
                 :type int :default 100 :help "Number of images in each mini-batch")
  (.add-argument parser "--epoch" "-e"
                 :type int :default 20 :help "Number of sweeps over the dataset to train")
  (.add-argument parser "--frequency" "-f"
                 :type int :default 1 :help "Frequency of taking a snapshot")
  (.add-argument parser "--gpu" "-g"
                 :type int :default -1 :help "GPU ID (negative value indicates CPU)")
  (.add-argument parser "--out" "-o"
                 :default "result" :help "Directory to output the result")
  (.add-argument parser "--resume" "-r"
                 :default "" :help "Resume the training from snapshot")
  (.add-argument parser "--unit" "-u"
                 :type int :default 1000 :help "Number of units")
  (setv args (.parse-args parser))

  (print (.format "GPU: {}" args.gpu))
  (print (.format "# unit: {}" args.unit))
  (print (.format "# Minibatch-size: {}" args.batchsize))
  (print (.format "# epoch: {}" args.epoch))
  (print "")

  ;; Set up a neural network to train
  ;; Classifier reports softmax cross entropy loss and accuracy at every
  ;; iteration, which will be used by the PrintReport extension below.

  (setv model (L.Classifier (MLP args.unit 10)))

  (when (>= args.gpu 0)
    (.use (chainer.cuda.get-device args.gpu)) ; Make a specified GPU current
    (.to-gpu model))                          ; Copy the model to the GPU
  
  ;; Setup an optimizer
  (setv optimizer (chainer.optimizers.Adam))
  (.setup optimizer model)

  ;; Load the MNIST dataset
  (setv mnist-dataset (chainer.datasets.get-mnist)
        train (first mnist-dataset)
        test  (second mnist-dataset)
        train-iter (chainer.iterators.SerialIterator train args.batchsize)
        test-iter  (chainer.iterators.SerialIterator test  args.batchsize
                                                     :repeat False :shuffle False))
  
  ;; Set up a trainer
  (setv updater (training.StandardUpdater train-iter optimizer :device args.gpu)
        trainer (training.Trainer updater [args.epoch "epoch"] :out args.out))

  ;; Evaluate the model with the test dataset for each epoch
  (trainer.extend (extensions.Evaluator test-iter model :device args.gpu))

  ;; Dump a computational graph from 'loss' variable at the first iteration
  ;; The "main" refers to the target link of the "main" optimizer.
  (trainer.extend (extensions.dump-graph "main/loss"))

  ;; Take a snapshot for each specified epoch
  (setv frequency (if (= args.frequency -1) args.epoch (max 1 args.frequency)))
  (trainer.extend (extensions.snapshot) :trigger [frequency "epoch"])
  
  ;; Write a log of evaluation statistics for each epoch
  (trainer.extend (extensions.LogReport))

  ;; Save two plot images to the result dir
  (when (extensions.PlotReport.available)
    (trainer.extend (extensions.PlotReport ["main/loss" "validation/main/loss"]
                                           "epoch" :file_name "loss.png"))
    (trainer.extend (extensions.PlotReport ["main/accuracy" "validation/main/accuracy'"]
                                           "epoch" :file_name "accuracy.png")))

  ;; Print selected entries of the log to stdout
  ;; Here "main" refers to the target link of the "main" optimizer again, and
  ;; "validation" refers to the default name of the Evaluator extension.
  ;; Entries other than 'epoch' are reported by the Classifier link, called by
  ;; either the updater or the evaluator.
  (trainer.extend (extensions.PrintReport
                   ["epoch" "main/loss" "validation/main/loss" "main/accuracy"
                            "validation/main/accuracy" "elapsed_time"]))

  ;; Print a progress bar to stdout
  (trainer.extend (extensions.ProgressBar))

  (if args.resume
    ;; Resume from a snapshot
    (chainer.serializers.load_npz args.resume trainer))
  
  ;; Run the training
  (trainer.run)
  )

(if (= --name-- "__main__") (main))

2017-05-11: ゼロから始めるHy(hylang)

Hy(hylang)について

HyはPythonのVM上で動くLisp方言で、Clojureによく似た構文を持つ。他のLispと同様にREPLでのインタラクティブな開発ができる。Lispなので本物のマクロが使える。

Pythonと高度な互換性がある。HyからPythonを呼ぶこともPythonからHyを呼ぶこともできる。

構文解析のオーバーヘッドが多少重いものの、Pythonのバイトコードへのコンパイラが付属しているので実用的には問題ない。またHyからPythonのソースコードへのトランスパイラも付属している。現状だとgensymの生成するシンボルをPythonのシンボルに変換できていないのでgensymを使うマクロがあるとエラーになってしまう。

Hyのインストール

Python2でも3でも動くようだが、Python3の方がサポートされている機能が多い。 最新版のHyではPython3.3以降のみがサポートされる。自分はpyenvで最新のPython3.6.0をインストールした。それから以下のようにするとGithub上のHyの最新版がインストールできる。

pip install git+https://github.com/hylang/hy.git

Emacsのhy-mode

簡易的なものだが、Emacs用のhy-modeが用意されている。MELPAから入るので、.emacsにpackage-install用の設定を書いて、

(require 'package)
(add-to-list 'package-archives 
	     '("melpa" . "http://melpa.milkbox.net/packages/") t)
(package-initialize)

しかるのちにM-x package-install RET hy-modeすればインストールされる。

.emacsの設定例

;;; Hy-mode

(autoload 'hy-mode "hy-mode"
  "Mode for editing Hylang source files"
  t)

(setq auto-mode-alist (append '(("\\.hy$" . hy-mode)) auto-mode-alist))

(add-hook 'hy-mode-hook
	  (lambda ()
            (setq hy-font-lock-keywords
                  (append '(("(\\|)" . paren-face))
                          hy-font-lock-keywords))))

REPLの起動

.hy拡張子のファイルを開くと、hy-modeになる。この状態でM-x inferior-lispとすると、*inferior-lisp*バッファにREPLが表示される。

以降はC-x eで式ごとに評価したり、C-c lでファイルをロードしたりできるようになる。

hy-mode.png

Hyのチュートリアル

Hyのチュートリアルを見ながら色々試してみた。

関数定義

docstringの位置が違う以外はほぼClojureと同じ(Clojureは何故あの位置なのだろう・・・)。

(defn fact [n]
  "Docstring"
  (if (= n 0)
    1
    (* n (fact (- n 1)))))

引数オプション(ラムダリスト)

関数やマクロのラムダリストには&optional&rest&kwargsを指定できる。

&optionalはCommon Lispのオプショナル引数とキーワード引数が混ざったようなもので、キーワードを指定しないとオプショナル引数のように振る舞うが、キーワードの名前を直接指定することで最初のキーワード引数を飛ばして二番目のキーワードだけを指定するようなことができる。

(defn optional-arg [pos1 pos2 &optional keyword1 [keyword2 42]]
  [pos1 pos2 keyword1 keyword2])

(optional-arg 'pos1 'pos2 'key1)         ; => ['pos1', 'pos2', 'key1', 42]
(optional-arg 'pos1 'pos2 :keyword2 420) ; => ['pos1', 'pos2', None, 420]

&restは可変長の引数を取り、それを関数/マクロの中で1つのリスト(実際にはPythonのタプル型)として扱える。 例えば次のコードは可変長の引数を取ってその総和を取る関数の定義になる。

(defmacro incf [var &optional [diff 1]]
  `(setv ~var (+ ~var ~diff)))

(defn plus [&rest args]
  (let ((sum 0))
    (for [i args] (incf sum i))
    sum))

&kwargsは&restのキーワード版みたいなもので、&restがタプルに対応したのに対して、&kwargsはディクショナリ型に対応する。

(defn some-func [foo bar &rest args &kwargs kwargs]
  (import pprint)
  (pprint.pprint (, foo bar args kwargs)))

(some-func 'foo 'bar 'arg1 'arg2)
;; => ('foo', 'bar', ('arg1', 'arg2'), {})

(some-func 'foo 'bar 'arg1 'arg2 :key1 'val1 :key2 'val2)
;;=> ('foo', 'bar', ('arg1', 'arg2'), {'key1': 'val1', 'key2': 'val2'})

変数

変数の宣言と代入の区別はないように見える。それぞれにdefsetvが用意されているが多分同じもので、単なる代入演算子( Pythonの= )だと思われる。

(def x 10)
(setv y 20)
(+ x y) ; => 30

PythonはLisp-1なので、HyもLisp-1である。すなわち、関数と変数で名前空間が分かれていない。従って、無名関数を使ってSchemeっぽい関数定義もできる。無名関数はClojureと同様にfnで作れる。

(def fact2
  (fn [n]
    (if (= n 0)
      1
      (* n (fact2 (- n 1))))))

letが無い!? → マクロで定義する

どうもletが用意されていないように見える。Pythonとの対応的には全てをsetvでやれということなのだろうか。とはいえletはラムダ式のシンタックスシュガーでしかないことを思い起せば簡単に定義できる。(参考: On Lisp — 古典的なマクロ)

マクロ定義はCommon LispやClojureのdefmacroとほとんど同じになっている。バッククォートとアンクォートを使うところも同じだ。アンクォートにはClojureと同じく~、スプライシングアンクォート(リストの埋め込み)には~@を使う。

マクロの展開形を確認するためにはmacroexpandを使う。また、変数捕捉を避けるためにgensymで重複しないことが保証されているシンボルを作れる。

(defmacro let [var-pairs &rest body]
  (setv var-names (list (map first  var-pairs))
        var-vals  (list (map second var-pairs)))
  `((fn [~@var-names] ~@body) ~@var-vals))

(macroexpand '(let ((one 1)
                    (two 2)
                    (three 3))
                (print one)
                (print two)
                (print three)))
;; (('fn' ['one' 'two' 'three'] ('print' 'one') ('print' 'two') ('print' 'three')) 1 2 3)

(let ((one 1)
      (two 2)
      (three 3))
  (print one)
  (print two)
  (print three))
;; 1
;; 2
;; 3

制御構造

まず真理値だが、基本Pythonと同じで、NoneとFalse、数字の0、空のシーケンス、空のディクショナリは偽、それ以外は真になる。それにしても0が偽になるのはどうかと思う。

(if (or None False 0 0.0 '() [] {}) 't 'f) ; => 'f'

condはClojureと少し違って、角括弧が必要なことに注意。

(let ((i 20))
  (cond [(> i 30) (print "That variable is too big!")]
        [(< i 10) (print "That variable is too small!")]
        [True     (print "That variable is jussssst right!")]))
;; That variable is jussssst right!

副作用を目的とした順次実行にはClojureと同様にdoを使う。これはCommon Lispのprogn、Schemeのbeginに相当する。Common LispとSchemeではdoはループのためのマクロなのでややこしい。

ifでは評価される部分に式が一つしか書けないので、複数の処理を書くためにはdoを使って処理をまとめる必要がある。

(if True
  (do
    (print "this is if true")
    (print "and why not, let's keep talking about how true it is!"))
  (print "this one's still simply just false"))
;; this is if true
;; and why not, let's keep talking about how true it is!

上の例でifにelse部分が無い場合は、whenを使った方が意図が明確になる。(副作用を目的としていることが分かるので)

(when True
  (print "this is if true")
  (print "and why not, let's keep talking about how true it is!"))
;; this is if true
;; and why not, let's keep talking about how true it is!

繰り返しにはforを使う。これはPythonのforそのままで、局所変数とシーケンスの組に続いてループ本体を書く。

(for [i (range 4)]
  (print (+ "'i' is now at " (str i))))

(for [i '(0 1 2 3)]
  (print (+ "'i' is now at " (str i))))

(for [i [0 1 2 3]]
  (print (+ "'i' is now at " (str i))))
  
;; 'i' is now at 0
;; 'i' is now at 1
;; 'i' is now at 2

関数定義のところで再帰でfactを書いたが、(fact 1000)とかにすると再帰が深すぎるといってエラーになってしまう。Hyは末尾再帰最適化はしないが、Clojureと同様にloop/recurマクロがあるので末尾再帰呼び出しを単純なループに変換できる。

(require [hy.contrib.loop [loop]])

(defn fact3 [n]
  (loop [[cnt 1] [acc 1]]
        (if (= cnt n)
          acc
          (recur (inc cnt) (* acc cnt)))))

(fact3 1000) ; => 4023872600770937735437024339 ...

Pythonの呼び出し

パッケージを読み込むimportはPythonとほとんど同じで、実行時に評価される。複数のimportを一つにまとめられる。

import sys
import os.path
(import sys os.path)

次に、PythonとHyでのfromasの使い方の対応を見る。

from os.path import exists, isdir as is_dir, isfile as is_file
from sys import *
import numpy as np
(import [os.path [exists
                  isdir :as dir?
                  isfile :as file?]]
        [sys [*]]
        [numpy :as np])

ただし、マクロは実行時でなくコンパイル時に評価されるので、マクロを含むパッケージはimportの代わりにrequireを使って読み込む。

(require [hy.contrib.loop [*]])    ; loop/recurはここ
(require [hy.extra.anaphoric [*]]) ; アナフォリックマクロのパッケージ

次に、NumPyの配列を作って属性にアクセスしてみる。

(def arr (np.array [[1 2 3]
                    [4 5 6]
                    [7 8 9]]))

arr.ndim  ; => 2
arr.size  ; => 9
arr.dtype ; => dtype('int64')

メソッド呼び出しは次の2つの書き方が両方通る。

(arr.sum)  ; => 45
(.sum arr) ; => 45

ほかにも色々やってみる。

;; スカラー倍
(* arr 3)
;; array([[ 3,  6,  9],
;;        [12, 15, 18],
;;        [21, 24, 27]])

;; アダマール積(要素積)
(* arr arr)
;; array([[ 1,  4,  9],
;;        [16, 25, 36],
;;        [49, 64, 81]])

;; 行列積
(np.dot arr arr)
;; array([[ 30,  36,  42],
;;        [ 66,  81,  96],
;;        [102, 126, 150]])

;; 一様乱数で100x100行列を作って行列積を取る
(import [numpy.random :as rand])

(def bigarr1 (rand.rand 100 100))
(def bigarr2 (rand.rand 100 100))

(np.dot bigarr1 bigarr2)

;; array([[ 28.38096367,  28.63420504,  28.01482173, ...,  27.26330009,
;;          25.56717227,  27.39401733],
;;        [ 25.26386499,  23.78039895,  22.81641922, ...,  24.37012314,
;;          22.31017675,  22.20606049],
;;        [ 24.79624411,  23.11758526,  24.45533016, ...,  24.47093385,
;;          22.3951194 ,  24.02735416],
;;        ..., 
;;        [ 25.65465691,  25.7403632 ,  23.54518075, ...,  24.36247407,
;;          21.92434498,  23.04834359],
;;        [ 22.37135022,  21.32717967,  21.92101116, ...,  20.93922527,
;;          20.07961519,  20.54109093],
;;        [ 27.50945536,  25.99902791,  25.73058543, ...,  25.71283456,
;;          23.86456424,  25.27311888]])

まとめ

  • Hy = Python + S式 + マクロ
    • Clojureライクな構文
    • Pythonのライブラリがそのまま使える
  • 次はChainerのサンプルコードをHyに翻訳してみる、かも

2017-04-25: cl-random-forest: Common Lispで高性能なランダムフォレストを実装した

はじめに

前の記事でCLMLのランダムフォレストを試したのだが、計算速度が遅くてとても実用レベルとは言えなかったので一から書き直すことにした。また先月のShibuya.lispのLispmeetupでも発表してきた。何をやっているかはこの発表のスライドで大体説明しているのだが、実際の使い方はデモでしか説明していなかったのでこの記事で説明する。

単純なランダムフォレストの構築ではCLMLの実装よりも50~60倍程度速くなっている。上記スライドにもベンチマークがあるが、RのRangerやPythonのscikit-learnなどの実装と比べても速く、精度も良い。

実装にあたって参考にしたのは、PRMLの14章とかこのへんのスライド

ランダムフォレストは分類だけでなく回帰にも使えるのだが、今のところ実装してるのは分類だけ。

論文: Global Refinement of Random Forest

単に普通のランダムフォレストを実装するだけでは面白くないので、加えて次の論文のアルゴリズムを実装している。

この手法は、学習済みのランダムフォレスト全体の構造を使って元のデータを高次元かつ疎なデータに変換する。そしてその疎なデータを線形分類器で再学習することによってより速く、より良い精度を出すというものだ。言い換えると、ランダムフォレストを非線形な特徴抽出器として使って線形分類器の入力を作っている。

更に、学習した線形分類器のパラメータ(重み)を使って寄与度の小さい枝を削除する枝刈りをすることもできる。そうすることによって最終的な予測精度をほとんど落とさずにモデルのサイズを大幅に縮小することができる。

この手法の再学習と枝刈りはランダムフォレスト内の決定木全体の情報を使って行なわれる。そのため、単純に各決定木の予測を平均化していた従来のランダムフォレストに対して、より相互補完的な予測の統合ができる。枝刈りに関しても、従来は個々の決定木ごとの情報量基準に従って枝刈りしていたが、この手法では森全体の中での重要度に応じて枝刈りするため、より効果が高い。

また、従来のランダムフォレストではバギング比率や分割の試行回数を減らすと性能が悪化してしまっていたが、この手法ではほとんど悪化しない(上記論文のFigure 3)。

ren2015-figure3.png

バギングに使うデータサイズを小さくすることができ、分割の試行回数も少なくできるということは、ランダムフォレストの構築にかかる時間を大幅に短縮できるということだ。そのためランダムフォレストの構築後に線形分類器の学習を行なったとしても、全体としての計算時間はむしろ短縮になる。

線形分類器は元論文ではLIBLINEARの線形SVMを使っているのだが、今回は前に実装したcl-online-learningのAROWを使う。その影響か元論文よりも若干の精度向上が見られた。

インストール

まだQuicklispに登録されていないので、local-projects以下に置くか、Roswellからインストールする。

ros install masatoi/cl-random-forest

最新のcl-online-learningが必要なので、Quicklispを最新の状態にするか、以下のように明示的にgithubから最新版をインストールするかする。

ros install masatoi/cl-online-learning masatoi/cl-random-forest

ランダムフォレストには少なからずメモリを使うので、処理系の起動オプションで動的に確保できるメモリ空間の最大値を大きめに取っておく。例えばSBCLに最大4GBのメモリを当てるには次のようにする。

sbcl --dynamic-space-size=4096

roswellの場合は

ros -Q dynamic-space-size=4096 run

使い方

データの用意

データセットは教師信号を表わす一次元fixnum配列と、データ本体を表わす二次元double-float配列からなる。このデータ行列における行が一つのデータに対応する。例えば、以下の図のような2次元4クラスのデータであれば、データセットのフォーマットはこうなる。教師信号(クラスID)は0スタートの整数であることに注意する。

clrf-dataset.png

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

(defparameter *datamatrix*
  (make-array '(11 2)
              :element-type 'double-float
              :initial-contents '((-1.0d0 -2.0d0)
                                  (-2.0d0 -1.0d0)
                                  (1.0d0 -2.0d0)
                                  (3.0d0 -1.5d0)
                                  (-2.0d0 2.0d0)
                                  (-3.0d0 1.0d0)
                                  (-2.0d0 1.0d0)
                                  (3.0d0 2.0d0)
                                  (2.0d0 2.0d0)
                                  (1.0d0 2.0d0)
                                  (1.0d0 1.0d0))))

決定木を作る

まずデータセットから一本の決定木を作ってみる。

(defparameter *n-class* 4) ; クラス数
(defparameter *n-dim* 2)   ; 特徴次元数

;; 決定木の構築(学習)
;; キーワードオプションは全部デフォルト値
(defparameter *dtree* 
  (clrf:make-dtree *n-class* *n-dim* *datamatrix* *target*
                   :max-depth 5          ; 決定木の最大深さ
                   :min-region-samples 1 ; 領域内の最小サンプル数、これを下回ったら分岐をやめる
                   :n-trial 10))         ; 分岐の試行回数

;; 予測: *datamatrix*の最初のデータを予測する
(clrf:predict-dtree *dtree* *datamatrix* 0) ; => 0

;; テスト
(clrf:test-dtree *dtree* *datamatrix* *target*)
; Accuracy: 100.0%, Correct: 11, Total: 11
; => 100.0, 11, 11

決定木は単なる構造体で、make-dtree関数で構築できる(決定木を構築した時点で学習は終わっている)。次に、predict-dtree関数で訓練に使ったデータセットの最初の行を入力にして、予測されるクラス番号を返している。この例では0が返っているが、*target*の最初の要素も0なので正解していることになる。データセット全体に対するテストはtest-dtree関数で行う。これは正答率、正答数、テストデータの数を多値で返す。

ランダムフォレストを作る

次に決定木の集合体であるランダムフォレストを作ってみる。

;; ランダムフォレストの構築
(defparameter *forest*
  (clrf:make-forest *n-class* *n-dim* *datamatrix* *target*
                    :n-tree 10         ; 決定木の数
                    :bagging-ratio 1.0 ; 元のデータに対してバギングで使うデータの比率
                    :max-depth 5
                    :min-region-samples 1
                    :n-trial 10))

;; 予測、テスト
(clrf:predict-forest *forest* *datamatrix* 0) ; => 0
(clrf:test-forest *forest* *datamatrix* *target*)
; Accuracy: 100.0%, Correct: 11, Total: 11
; => 100.0, 11, 11

ランダムフォレストも構造体で、内部に決定木のリストを持っている。ランダムフォレストはmake-forest関数で構築できる。キーワードオプションは決定木とほぼ同じだが、決定木の数n-treeとバギング比率bagging-ratioを指定するところが異なる。

バギングとは元のデータセットからブートストラップサンプリング(重複を許すランダムサンプリング)して小さなデータセットを決定木ごとに作り、それを使って各決定木を学習することをいう。バギングによって決定木の多様性が増すことによって全体としての頑健性が増す。

決定木と同様に、予測・テストにはそれぞれpredict-forest関数、test-forest関数を使う。

Global refinement

ここまでは通常のランダムフォレストの話だったが、ここからランダムフォレスト全体の情報を使っての再学習の話になる。

まずは、構築したランダムフォレストを使って、元のデータセットを線形分類器への入力となるデータセットに変換する。具体的には、入力データが決定木のどの葉ノードに分類されるか調べ、実際に分類される葉を1、それ以外の葉を0とする。これを各決定木について行ない、それを一列に並べたものが変換後のデータとなる。

clrf-global-refinement.png

この新しいデータセットはランダムフォレスト中の全決定木の葉の総数が特徴数になる非常に高次元のデータだが、実際に値が入るのは各決定木につき1つだけであり、疎なデータになる。従って、各データについて値が入っているところのインデックスだけをテーブルに保持しておけばいい。これにはmake-refine-datasetを使う。

(defparameter *forest-refine-dataset* (clrf:make-refine-dataset *forest* *datamatrix*))

;; テーブルのサイズはデータ数×決定木の数
;; #(#(3 5 11 15 19 23 26 30 33 39)
;;   #(3 7 11 15 19 23 26 30 35 39)
;;   #(2 5  9 14 17 22 25 28 32 37)
;;   #(2 4  9 14 16 22 25 28 31 37)
;;   #(1 6 10 13 18 20 26 29 34 38)
;;   #(1 6 10 13 18 21 26 29 34 38)
;;   #(1 6 10 13 18 21 26 29 34 38)
;;   #(0 4  8 12 16 20 24 27 31 36)
;;   #(0 4  8 12 17 20 24 27 31 36)
;;   #(0 4  8 12 17 20 24 27 32 36)
;;   #(0 4  8 12 17 21 25 27 32 36))

次に、このデータセットから線形分類器を学習する。make-refine-learner関数で線形分類器のオブジェクトを作り、train-refine-learner関数で訓練し、test-refine-learner関数でテストする。この簡単な例では学習データとテストデータを分けていないが、実際にはテストデータも変換する必要があることに注意する。

データによっては学習が収束するまでにはtrain-refine-learnerを数回呼ぶ必要があるかもしれない。簡単な収束判定を実装したものがtrain-refine-learner-process関数で、テストデータに対する精度で収束を判定するので訓練データと一緒にテストデータも与える必要がある。

(defparameter *forest-refine-learner* (clrf:make-refine-learner *forest*))
(clrf:train-refine-learner *forest-refine-learner* *forest-refine-dataset* *target*)
(clrf:test-refine-learner  *forest-refine-learner* *forest-refine-dataset* *target*)
;; Accuracy: 100.0%, Correct: 11, Total: 11

;; 収束判定付き訓練
(clrf:train-refine-learner-process
 *forest-refine-learner*
 *forest-refine-dataset* *target*  ; 訓練データ
 *forest-refine-dataset* *target*) ; テストデータ

データ1つを予測したい場合は次のようにする。

 (clrf:predict-refine-learner *forest* *forest-refine-learner* *datamatrix* 0)

Global pruning

さて、Global refinementの学習結果を利用して、重要度の小さい葉ノードを刈り込むことができる(pruning)。例えば、あるランダムフォレストから全体の10%の葉ノードを削除したい場合は以下のようにする。

;; ランダムフォレストを破壊的に枝刈り
(clrf:pruning! *forest* *forest-refine-learner* 0.1)

;; 線形分類器を再学習
(setf *forest-refine-dataset* (clrf:make-refine-dataset *forest* *datamatrix*))
(setf *forest-refine-learner* (clrf:make-refine-learner *forest*))
(clrf:train-refine-learner *forest-refine-learner* *forest-refine-dataset* *target*)
(clrf:test-refine-learner  *forest-refine-learner* *forest-refine-dataset* *target*)

pruning!関数を呼ぶことで*forest*が破壊的に変更される。枝刈りした後は線形分類器の再学習が必要になる。このランダムフォレストの枝刈り→線形分類器の再学習を繰り返していくことで、モデルがコンパクトになっていく。

並列化

ランダムフォレストの中の決定木は互いに独立しているので、簡単に並列化することができ、その並列化粒度も大きい。ランダムフォレストの構築make-forest、Global refinement用データへの変換make-refine-datasetは並列化することによってかなりの高速化が期待できる。

またマルチクラス分類の場合、Global refinement時の線形分類器は二値分類器をクラス数だけ用意して学習することになる。この二値分類器のセットも独立した計算なので並列化することができる。train-refine-learnertest-refine-learnertrain-refine-learner-processは並列化が効く。

並列化を有効化するには、lparallelのカーネル数を指定してやる必要がある。例えば自分の環境は4コアCPUなので、次のようにする。

;; 並列化有効化
(setf lparallel:*kernel* (lparallel:make-kernel 4))
;; 無効化
(setf lparallel:*kernel* nil)

今のところSBCL以外だと並列化はちゃんと動かないかも。

MNISTの例

もうちょっと実際的な例として、MNISTでやってみた例がここにある。

やっていることは上記の単純なデータでやっていることとほぼ同じで、テストデータに対して98.3%くらい出る。

次の図はMNISTでGlobal pruningをやってみた結果で、枝刈り回数に対する精度と葉ノード数を表している。葉ノード数が当初の1/10くらいになっても精度はあまり変わらないことが分かる。

clrf-mnist-pruning.png

ベンチマーク

有名なランダムフォレストの実装と比較してみる。比較に使ったのはPythonのパッケージscikit-learnとRのパッケージranger。それぞれ本体部分はCythonとC++で実装されていて、並列化にも対応しており、速いとされている。

推奨設定での精度と所要時間

  • scikit-learn: n-tree=100、他はデフォルト
  • ranger: 全部デフォルト
  • cl-random-forest: n-tree=500, bagging-ratio=0.1, max-depth=15, n-trial=特徴数の平方根

cl-random-forestはランダムフォレストの構築とGlobal refinementまでを行なう(枝刈りはしない)。

データはLIBSVMのサイトにあるデータセットを使う。

  #Train #Test #Feature #Classes
MNIST 60000 10000 784 10
letter 15000 5000 16 26
covtype 348607 232405 54 7
usps 7291 2007 256 10

covtypeは訓練データとテストデータが分かれていないので、全体をシャッフルしてから上の表にあるデータ数で分割した。

結果は以下の表のようになる。ほとんどの場合でcl-random-forestの方が速く、精度も良いことが分かった。

  scikit-learn ranger cl-random-forest
MNIST 96.95%, 41.72sec 97.17%, 69.34sec 98.29%, 12.68sec
letter 96.38%, 2.569sec 96.42%, 1.828sec 97.32%, 3.497sec
covtype 94.89%, 263.7sec 83.95%, 139.0sec 96.01%, 103.9sec
usps 93,47%, 3.583sec 93.57%, 11.70sec 94.96%, 0.686sec

まとめ

  • Common Lispでランダムフォレストの実装cl-random-forestをつくった
  • 論文「Global Refinement of Random Forest」のアルゴリズムを実装した
  • 有名なランダムフォレストの実装と比較した。速度でも精度でも勝っている
  • TODO

2017-04-24: Hello World

はてなダイアリーからの移行

はてなははてなダイアリーからはてなブログへ移行することを推奨しているのだが、これまでは移行するメリットを見出せず放りっぱなしにしていた。 しかし最近Chromeではてなダイアリーにアクセスすると異常に遅いことに気付いた。どうもメンテナンスされていないJavascriptが原因っぽい。それでいい加減に他のブログシステムに移行した方がいいと思うようになった。

日々のメモはEmacsのorg-modeで取っているので、できればブログもorg-modeで編集できた方がいいのだが、pandocなどでorg-modeからmarkdownに変換する方法もあるということなので、メジャーな方法、Jekyll bootstrap+Github pagesでいくことにした。これだとレスポンスも速いし情報源もたくさんある。

まずやったことはrbenvのインストール、最新のrubyをインストール (ruby2以降でないとJekyllがインストールできないので)。それからgithubでブログサイトのレポジトリをつくった。基本的にはここ(Jekyll QuickStart)のサイト名の部分だけを差し替えて実行すればいい。

それから_config.ymlをいじってサイト情報を書き込んだりシンタックスハイライトを有効化したりする。

テーマの変更

色々見たが、最初から入っているTwitterテーマを使うことにした。 ~/masatoi.github.io/assets/themes/twitter/css/style.cssをいじると反映される。シンタックスハイライトのCSSスタイルはpygentsのテーマから適当なものを選んで叩き台にする。

なお、highlighterはrougeでないとコミットしたときにgithubからお叱りのメールが来る。

markdown-modeのインストール

Emacsのpackageからmarkdown-modeを入れる。これでmarkdownのシンタックスハイライトとC-c C-c pでEmacs規定のブラウザでプレビューができる。もっとも、Jekyll serveでファイルの変更を監視していればmdファイルを保存した段階でlocalhost:4000でプレビューが見られる。

コードのテスト

(defun split-arr3 (pred arr true-array false-array)
  (declare (optimize (speed 3) (safety 0))
           (type (simple-array fixnum) arr true-array false-array))
  (let ((true-index 0)
        (false-index 0))
    (declare (type fixnum true-index false-index))
    (loop for elem fixnum across arr do
      (cond ((funcall pred elem)
             (setf (aref true-array true-index) elem)
             (incf true-index))
            (t
             (setf (aref false-array false-index) elem)
             (incf false-index))))
    (values true-array true-index false-array false-index)))
    
(defun main ()
  (print "common lisp start.")
  (loop for count from 1 to *n*
        collect (sumup3 count))
  (print "common lisp end."))

defunloopも組み込みマクロと同列に表示されてしまうのがかなり納得いかないのだが、とりあえずそれらしくなったからよしとする。

コメントシステム

コメントにはDisqusを使えるが、なんか重くなるしそもそもはてなダイアリー時代にもコメントがついたことなんてほとんどないので無しでいこうと思う。はてなブックマークへの誘導を用意すれば一応コメントも残せるはず。

インラインHTMLのテスト

画像のテスト

alt text

画像に限らず、サイトのディレクトリ直下にファイルを置いておけばURLから普通に参照できるようだ。