ブログトップ 記事一覧 ログイン 無料ブログ開設

翡翠はコンピュータに卵を生むか

2017-02-13

CLMLのランダムフォレストを試してみる

ランダムフォレストは決定木ベースの分類/回帰モデルで、ニューラルネットやSVMなどと同様に非線形モデルなので線形分離不可能な問題にも使える。SVMはデータ数に対して指数的に計算時間がかかる一方、ランダムフォレストはデータ数をnとしてn*log(n)のオーダであり、軽い。また、SVMは基本的に二値分類器なので複数の学習器を組み合わせてマルチクラス分類することが多いが、ランダムフォレストは元からマルチクラスに対応していて追加のコストがかからない。さらに並列化も簡単にできるなど、扱いやすい性質をたくさん備えているので現実世界でよく使われている。

参考にしたランダムフォレストについての記事いろいろ

CLMLにランダムフォレストの実装があったので試してみる。

CLMLはCommon Lisp用の機械学習パッケージ集であり、Quicklispからインストールできる。ただし、あらかじめ処理系のdynamic-space-sizeを2560以上にして起動しておく必要があることに注意する。

先に必要ライブラリを読み込んでおく。

(ql:quickload :clml)
(ql:quickload :cl-online-learning)
(ql:quickload :alexandria)

データの読み込み

まずはデータの読み込みの例。サンプルデータをネットからダウンロードしてきて読み込む。

(defparameter *syobu*
  (clml.hjs.read-data:read-data-from-file 
   (clml.utility.data:fetch "https://mmaul.github.io/clml.data/sample/syobu.csv")
   :type :csv :csv-type-spec '(string integer integer integer integer)))

;; CL-USER> *syobu*
;; #<CLML.HJS.READ-DATA:UNSPECIALIZED-DATASET >
;; DIMENSIONS: 種類 | がく長 | がく幅 | 花びら長 | 花びら幅
;; TYPES:      UNKNOWN | UNKNOWN | UNKNOWN | UNKNOWN | UNKNOWN
;; NUMBER OF DIMENSIONS: 5
;; DATA POINTS: 150 POINTS

;; CL-USER> (aref (clml.hjs.read-data:dataset-points *syobu*) 0)
;; #("Setosa" 51 35 14 2)

ここから分かるように、データセットは特徴量の各次元の名前のリストと、データ点を表すベクタのベクタを持っている。例えば以下のようにして新たにデータセットを作れる。

(defparameter dataset1
  (clml.hjs.read-data:make-unspecialized-dataset
   '("class" "feat1" "feat2")
   (vector (vector "posi" 1 2)
           (vector "nega" -1 -2)
           (vector "posi" 10 20))))

データ点のベクタにラベルの文字列と数値が混在しているのが遅そう。

LIBSVMデータセットのデータを読み込む

今読んでるランダムフォレストのオンライン版の論文に出てくるものと同じデータセットで比較したいので、LIBSVMのデータセット集からmushroomsデータを読み込む。

cl-online-learningのread-data関数で読み込んでCLMLの形式に変換する。予測したい特徴の名前が後で必要になってくるので、クラスラベルに"class"と名前を付けて、残りの特徴名には単に連番を振っておく。

(defparameter *mushrooms-dim* 112)
(defparameter *mushrooms-train* (clol.utils:read-data "/home/wiz/datasets/mushrooms-train" *mushrooms-dim*))
(defparameter *mushrooms-test* (clol.utils:read-data "/home/wiz/datasets/mushrooms-test" *mushrooms-dim*))

(defun clol.datum->clml.datum (datum)
  (let ((label (if (> (car datum) 0) "posi" "nega")))
    (coerce (cons label (coerce (cdr datum) 'list)) 'vector)))

(defun clol.dataset->clml.dataset (dataset)
  (let ((datum-dim (length (cdar dataset))))
    (clml.hjs.read-data:make-unspecialized-dataset
     (cons "class" (mapcar #'(lambda (x) (format nil "~A" x)) (alexandria:iota datum-dim)))
     (map 'vector #'clol.datum->clml.datum dataset))))

(defparameter *mushrooms-train.clml* (clol.dataset->clml.dataset *mushrooms-train*))
(defparameter *mushrooms-test.clml* (clol.dataset->clml.dataset *mushrooms-test*))

ランダムフォレストの学習

lparallelを使って並列化しているので、まずカーネルサイズの設定をしておく。

(setf lparallel:*kernel* (lparallel:make-kernel 4))

上で作った *mushrooms-train.clml* データセットから学習する。予測対象のクラスラベルの特徴名をここで指定する。

(defparameter *forest* (clml.decision-tree.random-forest:make-random-forest *mushrooms-train.clml* "class"))
;; Evaluation took:
;;   81.367 seconds of real time
;;   269.523032 seconds of total run time (265.196590 user, 4.326442 system)
;;   [ Run times consist of 9.140 seconds GC time, and 260.384 seconds non-GC time. ]
;;   331.24% CPU
;;   276,010,481,282 processor cycles
;;   120,282,960,608 bytes consed

4コアCPUで計算しているが、かなり時間がかかってしまっている。

予測は、予測したいデータ点のラベル部分を "?" に置き換えたものを学習済みのモデルとともにpredict-forest関数に与える。例えば、テストセットの最初のデータ点を予測するには、

(aref (clml.hjs.read-data:dataset-points *mushrooms-test.clml*) 0)
;; #("nega" 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0 1.0d0 0.0d0 0.0d0 0.0d0 0.0d0
;;   0.0d0 1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0 0.0d0 0.0d0
;;   1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0 1.0d0 0.0d0 0.0d0 1.0d0
;;   1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0
;;   1.0d0 0.0d0 0.0d0 1.0d0 0.0d0 0.0d0 1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0
;;   0.0d0 0.0d0 1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0 0.0d0
;;   0.0d0 1.0d0 0.0d0 0.0d0 1.0d0 0.0d0 0.0d0 1.0d0 0.0d0 0.0d0 0.0d0 1.0d0 0.0d0
;;   0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0 0.0d0 0.0d0 0.0d0 0.0d0
;;   0.0d0 1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0 0.0d0 0.0d0)

(defparameter *query*
  #("?" 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0 1.0d0 0.0d0 0.0d0 0.0d0 0.0d0
    0.0d0 1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0 0.0d0 0.0d0
    1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0 1.0d0 0.0d0 0.0d0 1.0d0
    1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0
    1.0d0 0.0d0 0.0d0 1.0d0 0.0d0 0.0d0 1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0
    0.0d0 0.0d0 1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0 0.0d0
    0.0d0 1.0d0 0.0d0 0.0d0 1.0d0 0.0d0 0.0d0 1.0d0 0.0d0 0.0d0 0.0d0 1.0d0 0.0d0
    0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0 0.0d0 0.0d0 0.0d0 0.0d0
    0.0d0 1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0 0.0d0 0.0d0))

(clml.decision-tree.random-forest:predict-forest *query* *mushrooms-train.clml* *forest*)
; => "nega"

テストセットをまとめて予測するには、forest-validation関数を使う。

(clml.decision-tree.random-forest:forest-validation *mushrooms-test.clml* "class" *forest*)

;; Evaluation took:
;;   12.064 seconds of real time
;;   12.075773 seconds of total run time (12.044405 user, 0.031368 system)
;;   [ Run times consist of 0.337 seconds GC time, and 11.739 seconds non-GC time. ]
;;   100.10% CPU
;;   40,922,800,082 processor cycles
;;   10,874,892,336 bytes consed

;; ((("nega" . "posi") . 100) (("posi" . "posi") . 428) (("nega" . "nega") . 1596))
;; => Error rate: 0.04708098

この返り値は、例えば"nega"のものを"posi"と判定したケースが100件あったということを意味する。残りは正解しているので正答率は95.3%程度となる。論文では99%くらい出るとあるのでちょっと低い。なんで。Hivemailでも99%くらい出る模様。

メタパラメータの指定

メタパラメータとして指定できるものは、

  • 決定木の数 (デフォルトは500)
  • 元のデータセットから重複を許してサンプリングしたものをデータセットとして各決定木を学習するのだが、その際のサンプルサイズ (デフォルトは元のデータセット全体のサイズ)
  • 元の特徴量からサンプリングしたものを各決定木の特徴量として使うのだが、その際のサンプリングする特徴の数 (デフォルトはフルの特徴量を使うが、よく使われるのは元の特徴量のサイズの平方根)

これらを調整してより小さなランダムフォレストを作ると、

(defparameter *small-forest*
  (clml.decision-tree.random-forest:make-random-forest
   *mushrooms-train.clml* "class"
   :tree-number 100
   :data-bag-size (floor (/ (length *mushrooms-train*) 10))
   :feature-bag-size (floor (sqrt *mushrooms-dim*))))

;; Evaluation took:
;;   1.217 seconds of real time
;;   4.279771 seconds of total run time (4.191569 user, 0.088202 system)
;;   [ Run times consist of 0.141 seconds GC time, and 4.139 seconds non-GC time. ]
;;   351.68% CPU
;;   4,125,539,396 processor cycles
;;   2,282,142,096 bytes consed

;; => Error rate: 0.041902073

となってかなり高速になり、精度も若干上がった。

番外編: cl-online-learningでもやってみる

(defparameter arow-learner (clol:make-arow *mushrooms-dim* 0.1d0))
(loop repeat 10 do
  (clol:train arow-learner *mushrooms-train*)
  (clol:test arow-learner *mushrooms-test*))

;; Accuracy: 91.99623%, Correct: 1954, Total: 2124
;; Accuracy: 94.3032%, Correct: 2003, Total: 2124
;; Accuracy: 93.92655%, Correct: 1995, Total: 2124
;; Accuracy: 93.87947%, Correct: 1994, Total: 2124
;; Accuracy: 93.83239%, Correct: 1993, Total: 2124
;; Accuracy: 93.87947%, Correct: 1994, Total: 2124
;; Accuracy: 93.83239%, Correct: 1993, Total: 2124
;; Accuracy: 93.83239%, Correct: 1993, Total: 2124
;; Accuracy: 93.83239%, Correct: 1993, Total: 2124
;; Accuracy: 93.83239%, Correct: 1993, Total: 2124

;; Evaluation took:
;; 0.018 seconds of real time
;; 0.018714 seconds of total run time (0.014723 user, 0.003991 system)
;; 105.56% CPU
;; 64,026,632 processor cycles
;; 786,432 bytes consed

(defparameter scw-learner (clol:make-scw *mushrooms-dim* 0.999d0 0.001d0))
(loop repeat 10 do
  (clol:train scw-learner *mushrooms-train*)
  (clol:test scw-learner *mushrooms-test*))

;; Accuracy: 38.93597%, Correct: 827, Total: 2124
;; Accuracy: 66.00753%, Correct: 1402, Total: 2124
;; Accuracy: 69.96233%, Correct: 1486, Total: 2124
;; Accuracy: 81.40301%, Correct: 1729, Total: 2124
;; Accuracy: 88.46516%, Correct: 1879, Total: 2124
;; Accuracy: 93.36158%, Correct: 1983, Total: 2124
;; Accuracy: 96.65725%, Correct: 2053, Total: 2124
;; Accuracy: 97.834274%, Correct: 2078, Total: 2124
;; Accuracy: 98.44633%, Correct: 2091, Total: 2124
;; Accuracy: 98.06968%, Correct: 2083, Total: 2124

;; Evaluation took:
;;   0.038 seconds of real time
;;   0.037766 seconds of total run time (0.037766 user, 0.000000 system)
;;   100.00% CPU
;;   128,015,604 processor cycles
;;   2,064,384 bytes consed

このデータセットではSCW-Iが比較的良い性能を出した。過学習しがちなので早めに止める必要がある。

まとめ/感想

  • CLMLのランダムフォレストを試してみた
  • データセットが型指定しないunspecialized-datasetであり、遅い
  • あまり最適化されていないのでまだまだ速くなりそう
  • 精度も何故か良くない
  • CLMLの実装を叩き台にしてオンラインランダムフォレストを実装してみるつもり

2016-12-19

Common Lispが機械学習に向いていると考えるこれだけの理由

Lisp Advent Calendar 2016参加記事

ここ数年ディープラーニングの出現をきっかけにAIが再び盛り上がっているので、いよいよLispの復権があるかと思いきや、ないので(泣)、多少なりともLispに興味を持ってもらえるように、LispとAIの関係について私見を述べておこうと思う。Lispといっても色々あるが、この記事では主にCommon Lispの話になる。

Lispというとどうしても過去の記号処理的AIと結びつけられてしまい、機械学習を駆使するような現代のAIでは役に立たないように思われがちなのだが、これは大体誤解である。少なくともCommon Lispは現代的なAI開発に適した特徴を備えている。まず、AI実装のためのプログラミング言語に必要とされる特徴は何なのかを明らかにするために、AIの歴史から考えてみたい。

AIの歴史

初期の記号処理的AI(以降は記号AIと呼ぶ)ではLispやPrologが実装言語として広く採用されてきた歴史がある。記号処理とはその名の通り記号を操作対象とする処理のことで、具体的には、エキスパートシステム、数式処理、プログラミング言語のコンパイラなど、論理的な推論や構造の変換を伴うものが多い。それらの記号ベースのデータ構造はリストで表現されることが多かったので、リストの取り扱いを得意とするLispが採用されてきたわけだ。*1 記号AIは一定の成功をおさめたが、現実世界の複雑な問題に適用しようとするとフレーム問題や記号接地問題の壁にぶち当たった。結局、古典的な記号AIでは事物の特徴をどのように取り出し、抽象化し、関連付ければいいのか、それをどこまでやったらいいのかといったことを解決できなかったのだ。

現実世界でも適用可能な妥当なルールベースを作りこむことが難しいことが分かってきたので、その後は大量のデータからボトムアップ的に学習して自動的に知識表現を獲得してやろうという流れになり、ニューラルネットなどの機械学習手法が発達した。最近のディープラーニングなどはその延長線上にある。

機械学習は数値計算

ほとんどの機械学習のアルゴリズムは何らかの目的関数を最適化することによって学習を実現している。例えばニューラルネットの学習であれば、モデルからの予測値と実測値の誤差を目的関数として、その勾配(微分)を取り、誤差を最小化させる方向にモデルパラメータを少し動かすことを繰り返す。このことからも分かるように、記号AIのプログラミングはリスト処理が中心だったのに対して、機械学習AIのプログラミングは数値計算が中心となる。従って、複雑な数値計算をいかに簡潔に記述でき、いかに効率的に計算できるかがキモになる。

より高水準に。より低水準に。

コンピュータの性能向上にともない、最終的な実行速度よりも開発効率を重視するというトレンドが生まれた。研究者は短時間でアルゴリズムを実装して実験する必要があり、自然とAI分野ではRやMATLAB、Pythonなどの動的で高水準な言語に人気が集まった。これらは基本的には逐次実行のスクリプト言語であり、実行速度は遅い。

これに対する戦略として、これらの言語では外部の数値計算ライブラリを呼び出す何らかの方法を持っている。具体的には、計算処理をある程度まとまった単位にまとめて、Cなどの低水準言語で書かれた数値計算ライブラリにまるごと送り出す。例えば、ループでやるような処理をベクトルに対するマップ操作として記述するとか、ベクトルの計算ならば、小単位のバッチ処理にまとめて行列の演算に置き換えてから外部ライブラリに渡すなどのテクニックを用いることになる。

しかしいつでもこれができるとは限らない。例えばリアルタイムのオンライン処理には向かないし、そもそも逐次処理を必要とする学習アルゴリズムでは使えない。そういう場合は結局Cなどで書いてコンパイルしたものを呼び出す必要がある。*2 こうなってくると、実行速度を出すためにアルゴリズムの選択やプログラミングスタイルに制約がかかることになり、本末転倒ということになりかねない。

実はCommon Lispは数値計算に向いている

もともとLispの特徴として、柔軟で強力な抽象化機構を持ち、時としてプログラマ自身が構文すら自由に拡張できる高水準言語だということがあった。それに加えてCommon LispにはSBCLなどのオープンソースの高品質なネイティブコンパイラがあり、Cで書かれたライブラリに簡単にアクセスできる仕組みCFFIがある。さらに最近ではインターネットからライブラリをダウンロード、インストール、ロードまでを一貫してこなすライブラリ管理システムQuicklispや、処理系のバージョンを管理したりデプロイを行なえるツールRoswellなども出てきて、いよいよ便利になりつつある。

一部ではLispはスクリプト言語で遅くて数値計算には向かないというイメージがあるようだが、SBCLなどはJITコンパイラであるし、現在ではコンパイル型言語の中でもむしろ速い部類に入る。この話題の概要としては、Gauche作者のshiroさんのエントリがとても分かりやすかった。

速いコードを書くための具体的な方法論として参考になる情報源をいくつか挙げておくと、

Common Lispが速いのは、動的言語にも関わらず、最適化宣言と型宣言をオプションとして付けることができるためだ。動的言語が遅いと言われるのは、主にオブジェクトの型情報を実行時に判定する必要があるためで、極端な話、全ての変数に型宣言を付けて、最適化宣言で実行時の型チェックを切れば、Common LispコンパイラはCとほぼ同等のネイティブコードを吐き出す。さらに、これらの宣言はあらゆるレベルで付けることができるので、関数単位、ブロック単位で、本当に必要な部分にのみ集中して最適化することができる。

そして、往々にしてそのようなボトルネックはプログラム全体のごく一部にすぎない。大抵のCommon Lisp処理系には、プロファイラという、プログラムのどこの部分に計算時間が割かれているかを調べるツールが付いており、ボトルネックを探すのに役に立つ。また、コンパイルした関数を個別にディスアセンブルすることもできる。アセンブリコードを見ることによって、実行時の型変換がどこで起こっているかや、メモリアロケーションがどこで発生しているかなどが分かるので、これもチューニングの手掛かりになる。*3

以上をまとめると、Common Lispの開発サイクルは以下のようになる。

  1. 動的で高水準な言語機能を使って迅速にプロトタイピングする
  2. プロファイリングでボトルネックを探す
  3. 最適化宣言、型宣言を付けてコンパイル
  4. ディスアセンブルしてネイティブコードを見ながら最適化

これら全てが開発環境の中から出ることなく、対話的に行うことができる。これにより開発サイクルが短く、試行錯誤できる回数が多くなり、結果として最適なアルゴリズムを見つけ出すチャンスが増える。往々にして「Common LispでCよりも速いコードが書ける」というのは、最高速度のことというよりは、この開発サイクル自体の効率性が高いことを示している。

必要なら外部の数値計算ライブラリにもアクセスできる

Common Lispはそれ自体で速いとはいえ、巨大な行列同士の掛け算のように、マルチコアCPUやGPUの性能を極限まで引き出したいときには外部ライブラリに頼ることになる。CFFI経由で比較的簡単にCの共有ライブラリから関数を呼ぶことができるが、主な数値計算ライブラリにはラッパーライブラリが存在する。Intel MKLやOpenBLASのラッパーライブラリとしてはLLAがあり、CUDAによる行列計算のラッパーライブラリとしてはMGL-MATがある。

数値計算ライブラリに投げた後の処理はPythonのNumpyなどと同じなのでここでの速度差はあまりないが、FFIの部分でCommon Lispの方が若干速くなる。

素のCommon Lispの速さがあれば、機械学習で最も計算能力を必要とする学習時には外部の数値計算ライブラリを使い、予測時にはピュアCommon Lispのみを使うという選択肢も取れるようになる。そうすることでアプリケーションとして配布するときの外部の依存ライブラリを減らし移植性を増すことができる。

Common Lispの機械学習ライブラリ

Common Lispの機械学習ライブラリは数はあまり無いものの、重要なところは抑えているような気がする。とはいえまだまだ少ないので新規参入者求む!

まとめ

  • 機械学習AIに求められているのは迅速に開発できて、なおかつ実行速度も速いこと
  • Common Lispは基本的には高水準言語だが、ボトルネックになっているところを集中的に最適化することもできるので、低水準言語という側面も持つ
    • アルゴリズム選択の自由度が高く、ほとんどの処理をCommon Lisp内で完結できるので、実行速度と開発効率のバランスを取りやすい
    • 本当に重い処理は外部ライブラリに投げることもできる
  • 機械学習ライブラリも一応揃ってる
  • 新規参入しよう!

*1:Lispの由来はリスト処理(LISt Processing)から来ているという話がある

*2:一応、PythonであればCythonなどのネイティブコードにコンパイルする仕組みもある。

*3:例えば、新たにオブジェクトを生成して返すようなコードは、結果受取り用の変数を引数として与えてそれに代入させればメモリアロケーションは発生せず、ガベージコレクションの時間を削減できる。

2016-11-22

Common Lispによる線形分類器ライブラリcl-online-learningを書いた

去年、オンライン機械学習本(クマ本)を読んで線形分類器を実装する記事を書いたり、それらのアルゴリズムをまとめてcl-online-learningというライブラリを作ってLispmeetupで紹介したりした。

その後放置していたのだが、最近になってもはや使わないようなアルゴリズムは削除したり、疎ベクトルへの対応や、学習器のCLOSオブジェクトを単なる構造体にするなどの大きな変更をした。このあたりで一度ちゃんと紹介記事を書いておこうかと思う。

cl-online-learningの特徴は、

  • アルゴリズム: パーセプトロン、AROW、SCW-I (おすすめはAROW)
  • 二値分類、多値分類に対応 (one-vs-one、one-vs-rest)
  • データが密ベクトル、疎ベクトルのどちらの場合にも対応
  • 純Common LispでC/C++のライブラリ(AROW++)を上回る速度

インストール

local-projectsディレクトリにソースを展開する。

cd ~/quicklisp/local-projects/
git clone https://github.com/masatoi/cl-online-learning.git

あるいは、Roswellがインストールされているなら単に

ros install masatoi/cl-online-learning

データの読み込み

1つのデータはラベル(+1/-1)と入力ベクトルのペア(cons)で、データのシーケンスがデータセットとなる。 libsvm datasetsの形式のファイルからデータセットを作るには、read-data関数が使える。とりあえずデータはlibsvm datasetsの二値分類データからa1aを使うことにする。

(defpackage :clol-user
  (:use :cl :cl-online-learning :cl-online-learning.utils :cl-online-learning.vector))

(in-package :clol-user)

(defparameter a1a-dim 123)
(defparameter a1a-train (read-data "/path/to/a1a"   a1a-dim))
(defparameter a1a-test  (read-data "/path/to/a1a.t" a1a-dim))

(car a1a-train)
;; (-1.0d0
;;  . #(0.0d0 0.0d0 1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0 0.0d0
;;      0.0d0 1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0
;;      0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0
;;      0.0d0 0.0d0 1.0d0 0.0d0 0.0d0 1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0
;;      0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0
;;      0.0d0 0.0d0 0.0d0 1.0d0 0.0d0 0.0d0 1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0
;;      1.0d0 0.0d0 1.0d0 1.0d0 0.0d0 0.0d0 0.0d0 1.0d0 0.0d0 0.0d0 1.0d0 0.0d0
;;      0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0
;;      0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0
;;      0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0
;;      0.0d0 0.0d0 0.0d0))

モデル定義

学習器のモデルは単なる構造体で、make-系関数で生成できる。その際いずれもデータの次元数を必要とする。その他にAROWは1個、SCWは2個のメタパラメータを指定する必要がある。パーセプトロン、AROW、SCWのモデルをまとめて定義すると、

(defparameter perceptron-learner (make-perceptron a1a-dim))
(defparameter arow-learner (make-arow a1a-dim 10d0))        ; gamma > 0
(defparameter scw-learner  (make-scw  a1a-dim 0.9d0 0.1d0)) ; 0 < eta < 1 , C > 0

訓練

データ1個を学習するには各学習器のupdate関数を使う。AROWならarow-update関数。これにデータの入力ベクトルとラベルを与えることで、arow-learnerが破壊的に更新される。

(arow-update arow-learner (cdar a1a-train) (caar a1a-train))
;; #S(AROW
;;  :INPUT-DIMENSION 123
;;  :WEIGHT #(0.0d0 0.0d0 -0.04d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0) ...
;;  :BIAS -0.04d0
;;  :GAMMA 10.0d0
;;  :SIGMA #(1.0d0 1.0d0 0.96d0 1.0d0 1.0d0 1.0d0 1.0d0 1.0d0 1.0d0 1.0d0) ...
;;  :SIGMA0 0.96d0)

これをデータセット全体に対して行うのがtrain関数である。

(train arow-learner a1a-train)

予測

こうして学習したモデルを使って、ある入力ベクトルに対して予測を立てるには各学習器のpredict関数を使う。AROWならarow-predict関数。

(arow-predict arow-learner (cdar a1a-test))
;; 1.0d0

正解の値(caar a1a-test)が-1.0d0なのでここは外してしまっている。

これをテストデータ全体に対して行ない、正答率を返すのがtest関数である。

(test arow-learner a1a-test)
;; Accuracy: 84.44244%, Correct: 26140, Total: 30956

となって84%弱の精度が出ていることが分かる。

マルチクラス分類

データの読み込み (MNIST)

マルチクラス分類ではデータのラベルが+1/-1ではなく、0以上の整数になる。例えばlibsvm datasetsからMNISTのデータを落としてきて読み込んでみる。読み込みはread-data関数にmulticlass-pキーワードオプションをつけて呼び出す。

(defparameter mnist-dim 780)
(defparameter mnist-train (read-data "/home/wiz/tmp/mnist.scale" mnist-dim :multiclass-p t))
(defparameter mnist-test  (read-data "/home/wiz/tmp/mnist.scale.t" mnist-dim :multiclass-p t))
;; このデータセットはラベルが1からではなく0から始まるので1足しておく
(dolist (datum mnist-train) (incf (car datum)))
(dolist (datum mnist-test)  (incf (car datum)))

(car mnist-train)
;; (5 . #(0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 0.0d0 ...))
モデル定義

マルチクラス分類は二値分類器の組み合わせで実現する。組み合せ方には色々あるが、cl-online-learningではone-vs-oneとone-vs-restを用意している。一般にone-vs-oneの方が精度が高いが、クラス数の二乗に比例する二値分類器が必要になる。一方のone-vs-restはクラス数に比例する。

例えばone-vs-oneで、二値分類器としてAROWを用いる場合の定義はこうなる。

(defparameter mnist-arow
  (make-one-vs-one mnist-dim    ; データの次元
                   10           ; クラス数
                   'arow 10d0)) ; 二値分類器の型とそのパラメータ

この構造体に対しても二値分類のときと同じくone-vs-one-update、one-vs-one-predict関数でデータを一つずつ処理できるし、train、test関数でデータセットをまとめて処理できる。

訓練、予測

データセットを8周訓練する時間を計測し、テストを行うコードは以下のようになる。

(time (loop repeat 8 do (train mnist-arow mnist-train)))
(test mnist-arow mnist-test)
;; Evaluation took:
;;   3.946 seconds of real time
;;   3.956962 seconds of total run time (3.956962 user, 0.000000 system)
;;   100.28% CPU
;;   13,384,797,419 processor cycles
;;   337,643,712 bytes consed

;; Accuracy: 94.6%, Correct: 9460, Total: 10000
liblinearの場合

高速な線形分類器とされるliblinearで同じデータを学習してみる。

wiz@prime:~/tmp$ time liblinear-train -q mnist.scale mnist.model
real    2m26.804s
user    2m26.668s
sys     0m0.312s

wiz@prime:~/tmp$ liblinear-predict mnist.scale.t mnist.model mnist.out
Accuracy = 91.69% (9169/10000)

こちらはデータの読み込みなども含めた時間なのでフェアな比較ではないが、大まかにいってcl-online-learningの方が大幅に速いといえる。また精度もcl-online-learning(AROW + one-vs-one)の方が良い。ちなみにliblinearのマルチクラス分類はone-vs-restを使っているらしい。

疎なデータの分類

a1aのデータを見ると気付くのは、ほとんどの要素が0の疎(スパース)なデータであるということだ。例えば「単語が文書に出現する回数」のような特徴量は高次元かつスパースになる。これをそのまま扱うと空間計算量も時間計算量も膨れ上がってしまうので、このようなデータではデータの次元数の長さのベクタを用意するのではなく、非零値のインデックスと値のペアだけを保持しておけばいい。 cl-online-learning.vectorパッケージに定義されているsparse-vector構造体がそれで、インデックスのベクタと値のベクタをスロットに持つ。

(make-sparse-vector
 (make-array 3 :element-type 'fixnum :initial-contents '(3 5 10))
 (make-array 3 :element-type 'double-float :initial-contents '(10d0 20d0 30d0)))

;; #S(CL-ONLINE-LEARNING.VECTOR::SPARSE-VECTOR
;;    :LENGTH 3
;;    :INDEX-VECTOR #(3 5 10)
;;    :VALUE-VECTOR #(10.0d0 20.0d0 30.0d0))

疎ベクトルの形でデータセットを読み込むにはread-data関数にsparse-pキーワードオプションをつけて呼び出す。試しに、1355191次元という超高次元のデータセットnews20.binaryを読み込んでみる。1つのデータはラベルとsparse-vector構造体のペアになっていることが分かる。

(defparameter news20.binary-dim 1355191)
(defparameter news20.binary (read-data "/home/wiz/datasets/news20.binary" news20.binary-dim :sparse-p t))

(car news20.binary)
;; (-1.0d0
;;  . #S(CL-ONLINE-LEARNING.VECTOR::SPARSE-VECTOR
;;       :LENGTH 3645
;;       :INDEX-VECTOR #(0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
;;                       ...
;;                       3636 3637 3638 3639 3640 3641 3642 3643 3644)
;;       :VALUE-VECTOR #(0.01656300015747547d0 0.01656300015747547d0
;;                       ...
;;                       0.01656300015747547d0)))

データ中の非零値の数をヒストグラムにしてみるとこうなる。

(ql:quickload :clgplot)
(clgp:plot-histogram (mapcar (lambda (d) (clol.vector::sparse-vector-length (cdr d)))
                             news20.binary) 200 :x-range '(0 3000))

f:id:masatoi:20161122213909p:image

1355191次元といってもほとんどのデータが2000次元以下なので疎なデータであることが分かる。

これを学習するためには、二値分類器としてarowの代わりにsparse-arowを使う。同様にパーセプトロンやSCWにもスパース版がある。

(defparameter news20.binary.arow (make-sparse-arow news20.binary-dim 10d0))
(time (loop repeat 20 do (train news20.binary.arow news20.binary)))
(test news20.binary.arow news20.binary)
;; Evaluation took:
;;   1.588 seconds of real time
;;   1.588995 seconds of total run time (1.582495 user, 0.006500 system)
;;   [ Run times consist of 0.006 seconds GC time, and 1.583 seconds non-GC time. ]
;;   100.06% CPU
;;   5,386,830,659 processor cycles
;;   59,931,648 bytes consed

;; Accuracy: 99.74495%, Correct: 19945, Total: 19996
AROW++の場合

同じことをC++によるAROW実装のAROW++を使ってやってみる。

wiz@prime:~/datasets$ arow_learn -i 20 news20.binary news20.binary.model.arow 
Number of features: 1355191
Number of examples: 19996
Number of updates:  37643
Done!
Time: 9.0135 sec.

wiz@prime:~/datasets$ arow_test news20.binary news20.binary.model.arow 
Accuracy 99.915% (19979/19996)
(Answer, Predict): (t,p):9986 (t,n):9993 (f,p):4 (f,n):13
Done!
Time: 2.2762 sec.
liblinearの場合
wiz@prime:~/datasets$ time liblinear-train -q news20.binary news20.binary.model

real    0m2.800s
user    0m2.772s
sys     0m0.265s
wiz@prime:~/datasets$ liblinear-predict news20.binary news20.binary.model news20.binary.out
Accuracy = 99.875% (19971/19996)

なおAROW++もliblinearもベクトルの内部表現は疎ベクトルでやっている模様。

疎なデータの分類(マルチクラス)

マルチクラス分類の場合でも同じことができるので、MNISTでやってみる。MNISTも画像データではあるが、六割程度は0なので疎なデータといえる。

この場合は疎ベクトルかつマルチクラス分類なので、read-data関数にsparse-pとmulticlass-pの両方のオプションをつけてデータを読み込む。密なデータの時と同様にmake-one-vs-oneを使うが、その引数にsparse-arowを指定するところが異なる。あとは大体一緒。

(defparameter mnist-train.sp (read-data "/home/wiz/tmp/mnist.scale" mnist-dim :sparse-p t :multiclass-p t))
(defparameter mnist-test.sp  (read-data "/home/wiz/tmp/mnist.scale.t" mnist-dim :sparse-p t :multiclass-p t))
;; このデータセットはラベルが1からではなく0から始まるので1足しておく
(dolist (datum mnist-train.sp) (incf (car datum)))
(dolist (datum mnist-test.sp)  (incf (car datum)))

(defparameter mnist-arow.sp (make-one-vs-one mnist-dim 10 'sparse-arow 10d0))
(time (loop repeat 8 do (train mnist-arow.sp mnist-train.sp)))

;; Evaluation took:
;;   1.347 seconds of real time
;;   1.348425 seconds of total run time (1.325365 user, 0.023060 system)
;;   [ Run times consist of 0.012 seconds GC time, and 1.337 seconds non-GC time. ]
;;   100.07% CPU
;;   4,570,387,768 processor cycles
;;   337,618,400 bytes consed

となって、約3倍の高速化となっている。

まとめ

  • 純Common Lispで線形分類器を書いた
    • 学習器をCLOSオブジェクトではなく単なる構造体にしたり、ベクトル演算の実行時の型チェックを外す、学習器の構造体を破壊的に更新して一時的なデータ構造を作らないなどのチューニングにより訓練部分はかなり速い。
    • AROWとSCW-Iは共分散行列の対角成分だけを使う近似をしている
  • コマンドラインベースのliblinearなどとは立ち位置が違うかも
  • Common Lispにはmecab互換の形態素解析エンジンcl-igoもあるので、文書分類などに応用できるかも

2016-10-31

MGLによる2クラスロジスティック回帰

人工知能に関する断創録 - Kerasによる2クラスロジスティック回帰

音声でもとても参考にさせていただいたこちらのブログで今度からKerasの連載をするらしく、まず最初に最もシンプルな例としてロジスティック回帰の例が紹介されていたので、これをMGLでもやってみる。とはいえMGLに二値のクロスエントロピーがなかったのでマルチクラスのロジスティック回帰を2クラスに適用しているところがちょっと違う。

コード(Gist) mgl-logistic-regression.lisp

必要ライブラリの読み込み

mgl-userはMGLを使うためのユーティリティ集というか、サンプルコードを集めたもの。データの正規化などもここで定義してある。

データがCSVなのでfare-csvで読み込んでparse-numberでパースする。clgplotはgnuplotのフロント。mgl-userとclgplotはQuicklispから入らないのでlocal-projects以下に置く必要がある。

(ql:quickload :mgl-user) ; https://github.com/masatoi/mgl-user
(ql:quickload :fare-csv)
(ql:quickload :parse-number)
(ql:quickload :clgplot) ; https://github.com/masatoi/clgplot

データの読み込みと正規化

データ: https://raw.githubusercontent.com/sylvan5/PRML/master/ch4/ex2data1.txt

参考エントリと同様にPRMLのデータを読み込む。fare-csv:read-csv-fileで読み込むと文字列のリストになっているのでparse-numberで数値に直す。それからデータ点の構造体(datum)のベクタにする。

(in-package :mgl-user)

;;; データの読み込み
(defparameter data-list
  (mapcar (lambda (line)
            (mapcar #'parse-number:parse-number line))
          (fare-csv:read-csv-file "/home/wiz/tmp/ex2data1.txt")))

(defparameter dataset
  (let ((dataset (make-array (length data-list))))
    (loop for i from 0 to (1- (length data-list))
          for line in data-list
          do
       (setf (aref dataset i)
             (make-datum
                  :id i :label (nth 2 line)
                  :array (make-mat 2 :initial-contents (subseq line 0 2)))))
    dataset))

データの各次元の平均が0、標準偏差が1になるようにデータを正規化する。まずdatasetのコピーを作って、そのコピーを破壊的に正規化する関数dataset-normalize!を呼ぶ。

;;; データの正規化
(defparameter dataset-normal (copy-dataset dataset))
(dataset-normalize! dataset-normal)

データの可視化

正規化したデータをプロットする。

(defun plot-dataset (dataset)
  (let ((positive-data (remove-if-not (lambda (datum) (= (datum-label datum) 1)) dataset))
        (negative-data (remove-if-not (lambda (datum) (= (datum-label datum) 0)) dataset)))

    (clgp:plot-lists
     (list (loop for datum across positive-data collect (mref (datum-array datum) 1))
           (loop for datum across negative-data collect (mref (datum-array datum) 1)))
     :x-lists (list (loop for datum across positive-data collect (mref (datum-array datum) 0))
                    (loop for datum across negative-data collect (mref (datum-array datum) 0)))
     :style 'points)))

(plot-dataset dataset-normal)

f:id:masatoi:20161031221441p:image

ロジスティック回帰モデルの定義

隠れ層なしで入力層から直接出力層に接続するようなネットワークになる。活性化関数はシグモイド関数とする。

;;; モデル定義
(defparameter fnn-sigmoid
  (build-fnn (:class 'fnn :max-n-stripes 100)
    (inputs (->input :size 2))
    (f1 (->sigmoid inputs))
    (prediction (->softmax-xe-loss (->activation f1 :name 'prediction :size 2) :name 'prediction))))

出力層の活性化関数はソフトマックス関数で、出力次元が2次元になるところが元エントリとは異なる。

f:id:masatoi:20161031221440p:image

図にするとこんな感じで、2つのクラスへの所属確率というか、信頼度みたいなものを出力する。この信頼度は常に合計1になり、この値が高い方のクラスに分類されたと解釈する。

訓練

次のようにすると訓練が始まる。モデルが破壊的に変更され、学習済のモデルが返るが、表示は何もない。

(train-fnn-process fnn-sigmoid dataset-normal :n-epochs 100)

途中経過を出力したい場合はこのようにする。第三引数は本当ならテストデータを指定して、エポックが変わるタイミングで訓練データとテストデータでの正答率が表示される。

(train-fnn-process-with-monitor fnn-sigmoid dataset-normal dataset-normal :n-epochs 100)

出力はこのようになって、100エポック回すと最終的に91%くらいの正答率になる。

(中略)
2016-10-31 18:35:40: ---------------------------------------------------
2016-10-31 18:35:40: training at n-instances: 10000
2016-10-31 18:35:40: test at n-instances: 10000
2016-10-31 18:35:40: pred. train bpn PREDICTION acc.: 91.00% (10000)
2016-10-31 18:35:40: pred. train bpn PREDICTION xent: 3.456e-3 (10000)
2016-10-31 18:35:40: pred. test  bpn PREDICTION acc.: 91.00% (100)
2016-10-31 18:35:40: pred. test  bpn PREDICTION xent: 3.456e-3 (100)
2016-10-31 18:35:40: Foreign memory usage:
foreign arrays: 0 (used bytes: 0)
CUDA memory usage:
device arrays: 114 (used bytes: 400,112, pooled bytes: 0)
host arrays: 0 (used bytes: 0)
host->device copies: 202, device->host copies: 20,602

最適化はモーメンタムSGDでやっているが、他にADAMなども指定できる。その辺の指定はtrain-fnn-processなどの定義の中でやっているので、細かな調整をしたいときはこれらの定義をいじる必要がある。

重みを見てみる

clumpsなどでモデルの中身を覗くことができる。例えばこのようにしてバイアスと重みを表示してみる。

(let* ((f1-activation (aref (clumps fnn-sigmoid) 2))
       (bias (aref (clumps f1-activation) 0))
       (weight (aref (clumps f1-activation) 2)))
  (describe bias)
  (describe weight)
  (list bias weight))
#<->WEIGHT (:BIAS PREDICTION) :SIZE 2 1/1 :NORM 2.76652>
  [standard-object]

Slots with :INSTANCE allocation:
  NAME               = (:BIAS PREDICTION)
  SIZE               = 2
  NODES              = #<MAT 1x2 AF #2A((1.956225 -1.9562262))>
  DERIVATIVES        = #<MAT 1x2 A #2A((0.0 0.0))>
  DEFAULT-VALUE      = 0
  SHARED-WITH-CLUMP  = NIL
  DIMENSIONS         = (1 2)
#<->WEIGHT (F1 PREDICTION) :SIZE 4 1/1 :NORM 4.55551>
  [standard-object]

Slots with :INSTANCE allocation:
  NAME               = (F1 PREDICTION)
  SIZE               = 4
  NODES              = #<MAT 2x2 AF #2A((-2.4052894 2.4102905) (-2.1375985 2.1420684))>
  DERIVATIVES        = #<MAT 2x2 A #2A((0.0 0.0) (0.0 0.0))>
  DEFAULT-VALUE      = 0
  SHARED-WITH-CLUMP  = NIL
  DIMENSIONS         = (2 2)

(#<->WEIGHT (:BIAS PREDICTION) :SIZE 2 1/1 :NORM 2.76652>
 #<->WEIGHT (F1 PREDICTION) :SIZE 4 1/1 :NORM 4.55551>)

これらの中のNODESが上の図で示したユニット間の矢印に対応している。

予測をプロット

次のようにしてネットワークの出力を可視化してみる。

(defun plot-prediction (dataset fnn class)
  (let* ((min-x1 (loop for x across dataset minimize (mref (datum-array x) 0)))
         (max-x1 (loop for x across dataset maximize (mref (datum-array x) 0)))
         (min-x2 (loop for x across dataset minimize (mref (datum-array x) 1)))
         (max-x2 (loop for x across dataset maximize (mref (datum-array x) 1)))
         (x1-list (loop for x from min-x1 to max-x1 by 0.1 collect x))
         (x2-list (loop for x from min-x2 to max-x2 by 0.1 collect x)))
    (clgp:splot-list
     (lambda (x1 x2)
       (mref (predict-datum fnn
                            (make-datum :id 0 :label 0d0
                                        :array (make-mat 2 :initial-contents (list x1 x2))))
             class))
     x1-list x2-list :map t)))

(plot-prediction dataset-normal fnn-sigmoid 0)
(plot-prediction dataset-normal fnn-sigmoid 1)

クラス毎の信頼度を出してプロットしてみたのが次の図になる。

f:id:masatoi:20161031221439p:image

f:id:masatoi:20161031221437p:image

データのプロットと見比べてみると、決定境界のあたりに信頼度0.5が来ていることが分かる。

同じ方法でマルチクラスへ拡張できる。


MGLで回帰:多次元出力

シンプルな例で多次元出力ができるかテスト

(ql:quickload :mgl-user)
(ql:quickload :clgplot)

(in-package :mgl-user)

;;; 入力1次元、出力2次元のデータ
(defparameter sin-cos-data
  (let* ((data-size 10000)
         (data (make-array data-size)))
    (loop for i from 0 to (1- data-size) do
      (let* ((x (coerce (* pi (1- (random 2.0))) 'single-float))
             (y1 (cos x))
             (y2 (sin x)))
        (setf (aref data i)
              (make-regression-datum :id i
                                     :target (make-mat 2 :initial-contents (list y1 y2))
                                     :array  (make-mat 1 :initial-contents (list x))))))
    data))

;; モデル定義
(defparameter fnn
  (build-fnn (:class 'regression-fnn :max-n-stripes 100)
    ;; Input Layer
    (inputs (->input :size 1))
    (f1-activations (->activation inputs :name 'f1 :size 256))
    (f1 (->relu f1-activations))
    (f2-activations (->activation f1 :name 'f2 :size 256))
    (f2 (->relu f2-activations))
    (prediction-activations (->activation f2 :name 'prediction :size 2))
    ;; Output Lump: ->squared-difference
    (prediction (->loss (->squared-difference (activations-output prediction-activations)
                                              (->input :name 'targets :size 2))
                        :name 'prediction))))

;; 訓練
(train-regression-fnn-process-with-monitor fnn sin-cos-data :n-epochs 100)

;; 予測をプロット
(let* ((x-list (wiz:seq (- (- pi) 1) (+ pi 1) :by 0.1))
       (result
        (loop for x in x-list collect
          (let ((result-mat (predict-regression-datum
                             fnn
                             (make-regression-datum :id 0
                                                    :target (make-mat 2 :initial-element 1.0)
                                                    :array (make-mat 1 :initial-element x)))))
            (list (mref result-mat 0)
                  (mref result-mat 1))))))
  (clgp:plot-lists (list (mapcar #'car result)
                         (mapcar #'cadr result)
                         (mapcar #'sin x-list)
                         (mapcar #'cos x-list))
                   :x-lists (list x-list x-list x-list x-list)
                   :title-list '("prediction-1dim" "prediction-2dim" "sin(x)" "cos(x)")))

f:id:masatoi:20161031004847p:image

f:id:masatoi:20161031004846p:image

普通に学習できてる。

[-π,π]のデータで学習させているのでその範囲からはみ出すとズレる。

2016-06-13

Common Lispによる深層学習ライブラリMGL(2): 回帰問題の場合

前回はMNISTの分類問題を解いたが、今回は回帰問題の場合。

コード https://github.com/masatoi/mgl-user/blob/master/src/regression.lisp

データ点の表現とそれをまとめてモデルに設定する部分

分類と回帰で何が違うかというと、まず教師信号がラベルでなく実数ベクトルになるので、データ点の構造体を定義しなおす必要がある。それから分類では出力層のtargetにバッチサイズ分のリストを設定していたところを、回帰では教師信号を保持するlumpを別に作って、そのnodesスロットに教師信号ベクトルをバッチサイズ分まとめた行列を設定する。

;;; Data expression

(defstruct regression-datum
  (id nil :type fixnum)
  (target nil :type mat)
  (array nil :type mat))

;;;; Sampling, clamping, utilities

(defun sample-regression-datum-array (sample)
  (regression-datum-array (first sample)))

(defun clamp-regression-data (samples mat)
  (assert (= (length samples) (mat-dimension mat 0)))
  (map-concat #'copy! samples mat :key #'sample-regression-datum-array))

(defun sample-regression-datum-target (sample)
  (regression-datum-target (first sample)))

(defun clamp-regression-target (samples mat)
  (assert (= (length samples) (mat-dimension mat 0)))
  (map-concat #'copy! samples mat :key #'sample-regression-datum-target))

;;; Regression FNN
(defclass regression-fnn (fnn) ())

(defmethod set-input (samples (bpn regression-fnn))
  (let* ((inputs (find-clump 'inputs bpn))
         (targets (find-clump 'targets bpn)))
    (clamp-regression-data samples (nodes inputs))
    (clamp-regression-target samples (nodes targets))))

データの正規化

MNISTでは入力は[0,1]に正規化されていたので特に何もしていなかったが、データセットによっては入出力の各次元ごとに、平均を引いて標準偏差で割って平均0、分散1にする正規化をやった方がいいことがある。

;;; Copy dataset
(defun copy-regression-dataset (dataset)
  (let ((new-dataset (map 'vector (lambda (datum) (copy-regression-datum datum)) dataset)))
    (loop for new-datum across new-dataset
          for datum across dataset do
            (setf (regression-datum-array new-datum) (copy-mat (regression-datum-array datum))
                  (regression-datum-target new-datum) (copy-mat (regression-datum-target datum))))
    new-dataset))

;;; Normalize
(defun regression-dataset-average (dataset)
  (let* ((first-datum-array (regression-datum-array (aref dataset 0)))
         (first-datum-target (regression-datum-target (aref dataset 0)))
         (input-dim (mat-dimension first-datum-array 0))
         (output-dim (mat-dimension first-datum-target 0))
         (input-ave (make-mat input-dim))
         (output-ave (make-mat output-dim)))
    (loop for datum across dataset do
      (axpy! (/ 1.0 (length dataset)) (regression-datum-array datum) input-ave)
      (axpy! (/ 1.0 (length dataset)) (regression-datum-target datum) output-ave))
    (values input-ave output-ave)))

(defun regression-dataset-variance (dataset input-ave output-ave)
  (let* ((input-dim   (mat-dimension input-ave 0))
         (output-dim  (mat-dimension output-ave 0))
         (input-var   (make-mat input-dim))
         (output-var  (make-mat output-dim))
         (input-diff  (make-mat input-dim))
         (output-diff (make-mat output-dim)))
    (loop for datum across dataset do
      ;; input
      (copy! input-ave input-diff)
      (axpy! -1.0 (regression-datum-array datum) input-diff)
      (.square! input-diff)
      (axpy! (/ 1.0 (length dataset)) input-diff input-var)
      ;; output
      (copy! output-ave output-diff)
      (axpy! -1.0 (regression-datum-target datum) output-diff)
      (.square! output-diff)
      (axpy! (/ 1.0 (length dataset)) output-diff output-var))
    (values input-var output-var)))

(defun regression-dataset-normalize! (dataset &key test-dataset (noise-degree 1.0))
  (let* ((first-datum-array (regression-datum-array (aref dataset 0)))
         (input-dim (mat-dimension first-datum-array 0))
         (first-datum-target (regression-datum-target (aref dataset 0)))
         (output-dim (mat-dimension first-datum-target 0))
         (input-noise (make-mat input-dim :initial-element noise-degree))
         (output-noise (make-mat output-dim :initial-element noise-degree)))
    (multiple-value-bind (input-ave output-ave)
        (regression-dataset-average dataset)
      (multiple-value-bind (input-var output-var)
          (regression-dataset-variance dataset input-ave output-ave)
        (axpy! 1.0 input-noise input-var)
        (axpy! 1.0 output-noise output-var)
        (.sqrt! input-var)
        (.sqrt! output-var)
        (.inv! input-var)
        (.inv! output-var)
        (flet ((normalize! (datum)
                 (axpy! -1.0 input-ave (regression-datum-array datum))
                 (geem! 1.0 input-var (regression-datum-array datum) 0.0 (regression-datum-array datum))
                 (axpy! -1.0 output-ave (regression-datum-target datum))
                 (geem! 1.0 output-var (regression-datum-target datum) 0.0 (regression-datum-target datum))))
          (loop for datum across dataset do (normalize! datum))
          (if test-dataset
            (loop for datum across test-dataset do (normalize! datum)))
          'done)))))

ここでaxpy!が行列同士の足し算(というか alpha * X + Y の結果をYに代入する)、geem!が成分ごとの積を破壊的にやるMGL-MATの関数で、ほかにも成分ごとの平方根や逆数を代入する.sqrt!や.inv!など色々ある。これらの関数はCUDAによる高速化が効く。

学習部分とモニタリング

分類問題では評価は正答率によってだったが、回帰問題では予測したものと教師信号とのRMSE(平均二乗誤差の平方根)で測る。

learnerに持たせるmonitorには、測り方を意味するmeasurerとそれをまとめるcounterの2つのスロットがある。measurerに設定されているmgl-bp::costが予測と教師信号との差のL2ノルムを表す。counterには単純に平均を取るbasic-counterとRMSEで平均を取るrmse-counterがあり、今回はrmse-counterを設定する。

あとは大体分類と同じ。

;;; Monitoring
(defun log-regression-cost (optimizer learner)
  (when (zerop (n-instances optimizer))
    (report-optimization-parameters optimizer learner))
  (log-msg "train/test at n-instances: ~S (~A ephochs)~%" (n-instances optimizer)
           (/ (n-instances optimizer) (length (training optimizer))))
  (log-padded
   (let ((bpn (bpn learner))
         (monitors (monitors learner)))
     (append
      (monitor-bpn-results (make-sampler (training optimizer) :max-n 10000) bpn (list (car monitors)))
      (if (test optimizer)
        (monitor-bpn-results (make-sampler (test optimizer)) bpn (cdr monitors))))))
  (log-mat-room)
  (log-msg "---------------------------------------------------~%"))

(defun train-regression-fnn-with-monitor
    (fnn training &key test (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-regression-cost :period ,(length training))
                       (:fn reset-optimization-monitors
                            :period ,(length training)
                            :last-eval 0))))
         (measurer (lambda (instances bpn)
                     (declare (ignore instances))
                     (mgl-bp::cost bpn)))
         (monitors (cons (make-instance 'monitor
                            :measurer measurer
                            :counter (make-instance 'rmse-counter
                                        :prepend-attributes '(:event "rmse." :dataset "train")))
                         (if test
                           (list (make-instance 'monitor
                                    :measurer measurer
                                    :counter (make-instance 'rmse-counter
                                                :prepend-attributes '(:event "rmse." :dataset "test")))))))
         (learner (make-instance 'bp-learner :bpn fnn :monitors monitors))
         (dateset (make-sampler training :n-epochs n-epochs)))
    (minimize optimizer learner :dataset dateset)
    fnn))

(defun train-regression-fnn-process-with-monitor
    (fnn training &key test (n-epochs 30) (learning-rate 0.1) (momentum 0.9) without-initialize)
  (with-cuda* ()
    (repeatably ()
      (if (null without-initialize)
        (init-bpn-weights fnn :stddev 0.01))
      (train-regression-fnn-with-monitor
       fnn training :test test :n-epochs n-epochs :learning-rate learning-rate :momentum momentum)))
  (log-msg "End")
  fnn)

モデル定義

ネットワークの定義部分もほとんど分類問題と同じだが、出力層の定義の仕方が少し違う。回帰問題の場合、出力層の活性化関数は恒等写像なので、予測結果は最後の->activationオブジェクトの出力部分になり、activations-output関数でアクセスできる。

出力層では予測と教師信号の二乗誤差を表す->squared-differenceと->lossを組み合わせる。->squared-differenceの生成関数の引数には->activationオブジェクトの出力部分とtargetsという名前を持つ->inputオブジェクトを与える。この->inputオブジェクトにset-inputで教師信号を設定する。

例えば、入力2次元、出力1次元のネットワークならこのように書ける。

;;;; Activation access utilities
(defun activations-output (activations)
  (aref (clumps activations) 3))

(defun find-last-activation (bpn)
  (let ((clumps-vec (clumps bpn)))
    (loop for i from (1- (length clumps-vec)) downto 0 do
      (let ((clump (aref clumps-vec i)))
      (typecase clump
        (->activation (return clump)))))))

(defparameter fnn-regression
  (build-fnn (:class 'regression-fnn :max-n-stripes 100)
    ;; Input Layer
    (inputs (->input :size 2))
    (f1-activations (->activation inputs :name 'f1 :size 1200))
    (f1 (->relu f1-activations))
    (f2-activations (->activation f1 :name 'f2 :size 1200))
    (f2 (->relu f2-activations))
    (f3-activations (->activation f2 :name 'f3 :size 1200))
    (f3 (->relu f3-activations))
    (prediction-activations (->activation f3 :name 'prediction :size 1))
    ;; Output Lump: ->squared-difference
    (prediction (->loss (->squared-difference (activations-output prediction-activations)
                                              (->input :name 'targets :size 1))
                        :name 'prediction))))

予測

出力層の予測部分の与え方が違うだけでほとんど同じ

(defun predict-regression-datum (fnn regression-datum)
  (let* ((a (regression-datum-array regression-datum))
         (len (mat-dimension a 0))
         (input-nodes (nodes (find-clump 'inputs fnn)))
         (output-nodes (nodes (activations-output (find-last-activation 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))))

例: Rastrigin関数の近似

前に2個体分散GAの記事を書いたときに回帰の例として使ったRastrigin関数を近似してみる。

まずランダムな入力に対するRastrigin関数の値で訓練データセットを作り、格子状の入力でテストデータセットを作る。

それからデータを正規化する。訓練データの正規化パラメータと同じものでテストデータも正規化する必要がある。

最後に先に定義したfnn-regressionを使って訓練を実行する。

;;; Rastrigin function

(defun rastrigin (x-list)
  (let ((n (length x-list)))
    (+ (* 10 n)
       (loop for xi in x-list summing
	 (- (* xi xi) (* 10 (cos (* 2 pi xi))))))))

(defparameter *rastrigin-dataset*
  (let* ((n 10000)
         (arr (make-array n)))
    (loop for i from 0 to (1- n) do
      (let ((x (- (random 10.24) 5.12))
            (y (- (random 10.24) 5.12)))
        (setf (aref arr i) (make-regression-datum
                            :id (1+ i)
                            :target (make-mat 1 :initial-element (rastrigin (list x y)))
                            :array  (make-mat 2 :initial-contents (list x y))))))
    arr))

(defparameter *rastrigin-test*
  (let* ((n (* 64 64)) ; separate to 64x64 cells (by 0.16)
         (arr (make-array n)))
    (loop for i from 0 to 63 do
      (loop for j from 0 to 63 do
        (let ((x (- (* i 0.16) 5.04))
              (y (- (* j 0.16) 5.04)))
          (setf (aref arr (+ (* i 64) j))
                (make-regression-datum
                 :id (1+ (+ (* i 64) j))
                 :target (make-mat 1 :initial-element (rastrigin (list x y)))
                 :array  (make-mat 2 :initial-contents (list x y)))))))
    arr))

(defparameter *rastrigin-dataset-normal* (copy-regression-dataset *rastrigin-dataset*))
(defparameter *rastrigin-test-normal* (copy-regression-dataset *rastrigin-test*))
(regression-dataset-normalize! *rastrigin-dataset-normal*
                               :test-dataset *rastrigin-test-normal*
                               :noise-degree 0.0)

;; Run
(train-regression-fnn-process-with-monitor
 fnn-regression *rastrigin-dataset-normal*
 :test *rastrigin-test-normal* :n-epochs 500)

学習の進行過程の図。縦軸がRMSEで小さいほど良い。横軸は訓練データを何周したか。

f:id:masatoi:20160613024335p:image

左がオリジナルのテストデータで右が予測したもの。

f:id:masatoi:20160613024557p:image

ほぼ再現できている。

TODO

次はLSTMを使った時系列データの認識をやってみる