書籍「オンライン機械学習」を買ったのでCommon Lispで実装してみた。

オンライン機械学習 (機械学習プロフェッショナルシリーズ)
オンライン機械学習を買ったので、書いてあるアルゴリズムから線形識別器をいくつか試してみた。オンライン学習とは、データを一括処理するバッチ学習に対して、個々のデータを逐次処理する学習手法のこと。オンライン学習だと以下のようなうれしいことがある。

  • 実装が簡単: 基本的に実装は1ステップ毎の処理を実装するだけでほとんど終わる(スッキリ!)
  • リアルタイム処理に組込める: その時点までの観測を元にした学習結果が常に予測に使えるので、学習しながら予測もしなければならない応用に向いている
  • 収束が速く、データの性質が徐々に変化しても追従できる: これはリアルタイム処理に対して使うにはいい性質だが、偏ったデータやノイズに弱いという諸刃の剣でもある(そこで最適化対象に正則化項を入れるなどして悪影響を減らそうとする)。

データを一つ一つ処理するので、オンライン学習の計算量はデータ数に対して線形に増加する。バッチ学習のカーネル法などでは、データ数について2~3乗のオーダで計算量が増えていったりして、数万を超える大規模データには使えない。例えば非線形カーネルを使うSVMではかなり速いとされるSMOアルゴリズムでも時間計算量がO(mn^2)だが、線形カーネルに限ればSVMはオンラインに学習することができてO(mn)で済む (m: 入力次元、n: データ数)。
この本で扱っているほとんどの学習手法は線形分離可能を前提とするが、データを高次元化させたりして実用レベルの精度を出せることも多い。例えばTwitterのストリーミングAPIからデータを拾ってきて、単語のあるなしに対応する特徴量を持つ数万次元の疎なベクトルを入力として学習し、リアルタイムに何らかの予測を立てるような応用に適している。この本の後半にはそういった応用に向けて、疎なベクトル向けのベクトル演算の実装の仕方のようなものも解説されている。全体的に実装を指向した本という印象である。

実装言語はCommon Lisp

最近、機械学習界隈で話題に上るのはPythonなどのスクリプト言語が多いが、実際の学習アルゴリズム部分はCなどの高速な言語で実装されたライブラリ中にあって、まとまったデータをライブラリに渡しているだけの場合が多い気がする。まとまったデータを渡すバッチ処理ならそれでもいいかもしれないが、オンライン処理ではこの外部呼び出しのオーバーヘッドが問題となるし、そもそも単純なループ処理が遅いスクリプト言語にはあまり向かないように思える。
上述のようにオンラインで学習しながらその都度予測して何かをやるというアプリケーションを書くには、数値計算のような低水準なものからドメイン特化言語の定義のような高水準なものまでカバーする言語で書きたいのが人情である。というわけでCommon Lispを用いる。

コード

quicklispですぐ読み込める形にしたものをgithubに置いておく
インストールするには、まずシェルから

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

次にLisp処理系からquicklispでロードする。

(ql:quickload :cl-online-learning)

ベクトル演算は副作用あり

ベクトル演算はとりあえず密ベクトルに絞って考えて、CLML(Common Lisp Machine Learning)に付いてたユーティリティの一部を流用する。ベクトルはdouble-float型のsimple-arrayで表現する。ベクトル同士の足し算やスカラー倍などの操作において、結果受け取り用のベクトルを渡して、それを破壊的に変更することでベクトル演算の度にmake-arrayさせないようになっている。

(let ((v1 (make-array 3 :element-type 'double-float :initial-contents '(1d0 2d0 3d0)))
      (v2 (make-array 3 :element-type 'double-float :initial-contents '(10d0 20d0 30d0)))
      (result (make-array 3 :element-type 'double-float :initial-element 0d0)))
  (print (v+ v1 v2 result))
  (print (v-scale v1 3d0 result))
  (print (v+ result v2 result))
  (print result))

; #(11.0d0 22.0d0 33.0d0) 
; #(3.0d0 6.0d0 9.0d0) 
; #(13.0d0 26.0d0 39.0d0) 
; #(13.0d0 26.0d0 39.0d0) 

ベクトル演算の返り値は結果受け取り用ベクトルと同一オブジェクトなので、返り値を束縛した後に同じ結果受け取り用ベクトルを使って別のベクトル演算をやると値が変わってハマることになる(ハマった)。

予測とテスト部分

線形識別器なので、重みと入力ベクトルの内積を取って、それが正か負かでクラスを分ける。

;;; 符号関数
(defun sign (x)
  (if (> x 0d0) 1d0 -1d0))

;;; 線形識別器の決定境界
(defun f (input weight bias)
  (+ (inner-product weight input) bias))

;;; 線形識別器の予測
(defun predict (input weight bias)
  (sign (f input weight bias)))

どうせシーケンシャルアクセスしかしないので、データセットは教師信号と入力ベクトルのドット対のリストということにする。するとテストはこう書ける。

(defun test (test-data weight bias)
  (let ((len (length test-data))
	(n-correct (count-if (lambda (datum)
			       (= (predict (cdr datum) weight bias) (car datum)))
			     test-data)))
    (format t "Accuracy: ~f%, Correct: ~A, Total: ~A~%" (* (/ n-correct len) 100.0) n-correct len)))

訓練部分

とりあえず最も基本となるパーセプトロンと、線形SVMを確率的勾配法で更新する場合を実装してみる。
バイアス(決定境界fの切片)は、入力ベクトルの次元を1つ上げて(常に値が定数となる要素を入力ベクトルに追加する)、普通に重みを学習すると、重みベクトルの余計に1次元増えた要素がバイアスになっている。だが、この本によるとバイアスの更新は他の重みとは分けた方がいいとある。バイアスだけ学習率を変えるなどのヒューリスティクスがあるらしい。

;;; 3.3 パーセプトロン (アルゴリズム3.1)
;; 次のステップのweightとbiasを返す。なおweightは破壊的に変更される。
(defun train-perceptron-1step (input weight bias training-label)
  (if (<= (* training-label (f input weight bias)) 0d0)
    (if (> training-label 0d0)
      (values (v+ weight input weight) (+ bias 1d0))
      (values (v- weight input weight) (- bias 1d0)))
    (values weight bias)))

(defun train-perceptron-all (training-data weight bias)
  (loop for datum in training-data do
    (setf bias (nth-value 1 (train-perceptron-1step (cdr datum) weight bias (car datum)))))
  (values weight bias))

(defun train-perceptron (training-data)
  (let ((weight (make-dvec (length (cdar training-data)) 0d0))
	(bias 0d0))
    (train-perceptron-all training-data weight bias)))

;;; 3.6 サポートベクトルマシン (アルゴリズム3.3) 線形SVM + 確率的勾配法(SGD)
(defun train-svm-sgd-1step (input weight bias learning-rate regularization-parameter
			    training-label v-scale-result)
  (let* ((update-p (<= (* training-label (f input weight bias)) 1d0))
	 (tmp-weight
	  (if update-p
	    (v+ weight (v-scale input (* learning-rate training-label) v-scale-result) weight)
	    weight))
	 (tmp-bias (if update-p (+ bias (* learning-rate training-label)) bias)))
    (values
     (v-scale tmp-weight (- 1d0 (* 2d0 learning-rate regularization-parameter)) weight)
     (* tmp-bias (- 1d0 (* 2d0 learning-rate regularization-parameter))))))

(defun train-svm-sgd-all (training-data weight bias learning-rate regularization-parameter v-scale-result)
  (loop for datum in training-data do
    (setf bias
	  (nth-value 1 (train-svm-sgd-1step (cdr datum) weight bias
					    learning-rate regularization-parameter
					    (car datum) v-scale-result))))
  (values weight bias))

(defun train-svm-sgd (training-data learning-rate regularization-parameter)
  (let ((weight (make-dvec (length (cdar training-data)) 0d0))
	(bias 0d0)
	(v-scale-result (make-dvec (length (cdar training-data)) 0d0)))
    (train-svm-sgd-all training-data weight bias learning-rate regularization-parameter v-scale-result)))

データを用意して実験

ここまでで訓練、予測、テストができたので、libsvmのサイトからデータセットを取ってきて試してみる。

$ wget http://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/binary/a1a
$ wget http://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/binary/a1a.t
$ wget http://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/binary/a9a
$ wget http://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/binary/a9a.t

a1aは123次元、訓練データ1605個、a9aは123次元、訓練データ32,561個で、a9aの方はSVMで訓練させようとすると結構時間がかかる。これを読み込んで、教師信号と入力ベクトルのドット対のリストにする関数read-libsvm-dataを定義する。

(defun read-libsvm-data (data-path data-dimension)
  (let ((data-list nil))
    (with-open-file (f data-path :direction :input)
      (labels ((read-loop (data-list)
		 (let ((read-data (read-line f nil nil)))
		   (if (null read-data)
		     (nreverse data-list)
		     (let* ((dv (make-array data-dimension :element-type 'double-float :initial-element 0d0))
			    (d (ppcre:split "\\s+" read-data))
			    (index-num-alist
			     (mapcar (lambda (index-num-pair-str)
				       (let ((index-num-pair (ppcre:split #\: index-num-pair-str)))
					 (list (parse-integer (car index-num-pair))
					       (coerce (parse-number:parse-number (cadr index-num-pair)) 'double-float))))
				     (cdr d)))
			    (training-label (coerce (parse-integer (car d)) 'double-float)))
		       (dolist (index-num-pair index-num-alist)
			 (setf (aref dv (1- (car index-num-pair))) (cadr index-num-pair)))
		       (read-loop (cons (cons training-label dv) data-list)))))))
	(read-loop data-list)))))

(defparameter a1a-train (read-libsvm-data "/path//to/a1a"   123))
(defparameter a1a-test  (read-libsvm-data "/path//to/a1a.t" 123))
(defparameter a9a-train (read-libsvm-data "/path//to/a9a"   123))
(defparameter a9a-test  (read-libsvm-data "/path//to/a9a.t" 123))
;;;;; a1aに対する訓練とテスト

;; パーセプトロン
(multiple-value-bind (weight bias)
    (train-perceptron a1a-train)
  (test a1a-test weight bias))
; Accuracy: 81.806435%, Correct: 25324, Total: 30956

;; 線形SVM+SGD
(let ((learning-rate 0.01d0)
      (regularization-parameter 0.01d0))
  (multiple-value-bind (weight bias)
      (train-svm-sgd a1a-train learning-rate regularization-parameter)
    (test a1a-test weight bias)))
; Accuracy: 83.09859%, Correct: 25724, Total: 30956
;;;;; a9aに対する訓練とテスト

;; パーセプトロン
(multiple-value-bind (weight bias)
    (train-perceptron a9a-train)
  (test a9a-test weight bias))
; Accuracy: 79.988945%, Correct: 13023, Total: 16281

;; 線形SVM+SGD
(let ((learning-rate 0.01d0)
      (regularization-parameter 0.001d0))
  (multiple-value-bind (weight bias)
      (train-svm-sgd a9a-train learning-rate regularization-parameter)
    (test a9a-test weight bias)))
; Accuracy: 84.601685%, Correct: 13774, Total: 16281

となって、このデータでは線形SVMも通常のSVMとほぼ同等の性能を出すことが分かった。
訓練時間を比較してみると、

;;; パーセプトロン

(time (train-perceptron a1a-train))
; Evaluation took:
;   0.002 seconds of real time
;   0.002046 seconds of total run time (0.002046 user, 0.000000 system)
;   100.00% CPU
;   4,414,171 processor cycles
;   131,072 bytes consed

(time (train-perceptron a9a-train))
; Evaluation took:
;   0.011 seconds of real time
;   0.011445 seconds of total run time (0.011445 user, 0.000000 system)
;   100.00% CPU
;   26,044,485 processor cycles
;   2,850,816 bytes consed

;;; 線形SVM

(time (let ((learning-rate 0.01d0)
	    (regularization-parameter 0.01d0))
	(train-svm-sgd a1a-train learning-rate regularization-parameter)))
; Evaluation took:
;   0.004 seconds of real time
;   0.004468 seconds of total run time (0.004468 user, 0.000000 system)
;   100.00% CPU
;   9,648,051 processor cycles
;   393,216 bytes consed

(time (let ((learning-rate 0.01d0)
	    (regularization-parameter 0.001d0))
	(train-svm-sgd a9a-train learning-rate regularization-parameter)))
; Evaluation took:
;   0.023 seconds of real time
;   0.023348 seconds of total run time (0.023348 user, 0.000000 system)
;   100.00% CPU
;   53,203,526 processor cycles
;   8,192,000 bytes consed

となって、データ数に比例して計算時間が増えているのが分かる。SVMでは20秒くらいかかっていたのが一瞬で終わるのがすごい。これなら数GBクラスのデータセットでも扱えそうである。
今回はデータセットをリストで表現しているが、これが数GBクラスのサイズになってくるとメモリ上にとっておくのが大変になってくる。しかしオンライン学習ならデータをファイルから一つ一つ読み込みながら学習することもできるので、メモリ量による制約もほとんどない。低スペックな普通のPCでも本格的なビッグデータを解析できるという利点は大きい。

さらに高度な手法(AROW、SCW、ASGD)も今度暇なときにでも実装してみようと思う。

つづき 書籍「オンライン機械学習」を買ったのでCommon Lispで実装してみた。(AROW編)