Common Lispで行列演算: LLA(Lisp Linear Algebra)を使う
Common LispでCPUを使った行列演算を高速にやりたいときは、BLASやLAPACKといった、専用に書かれた外部ライブラリを呼ぶことになる。そのためのラッパーライブラリとして、LLA (Lisp Linear Algebra)がある。
LLAでは、普通のReference BLAS以外にもOpenBLAS(旧GotoBLAS)、ATLAS、MKLを使える。このうちマルチスレッド対応なのはATLASとOpenBLAS、MKLだが、ATLASはソースコードからビルドしないとオートチューンが効かない。Ubuntuでapt-getで入るATLASはシングルスレッドでしか動かないらしい（ubuntu 14.04でBLASを使う）。とはいえOpenBLAS と ATLAS の性能を R 上で比較するを見る限り、オートチューンできたとしてもATLASよりはOpenBLASを使った方が多分速い。
OpenBLASのインストール
OpenBLASはUbuntuならapt-getですぐに入る。
sudo apt-get install libopenblas-base sudo apt-get install libopenblas-dev
MKLのインストール
MKL(Math Kernel Library)はintel謹製の数値計算ライブラリで、intel製のマルチコアCPUの性能をフルに引き出すことができるとされる。
無償版はintelのページで個人情報を登録するとダウンロードできる。1GBくらいある。ダウンロード画面で表示される製品登録キーがインストール時に必要なのでメモっておくことを忘れずに。
インストール自体は、GUIインストーラがあるのでrootで起動して指示に従っていけばいい。その場合インストール先は/opt/intel/以下になる。
LLAのインストールと設定
LLAはquicklispから簡単にインストールできる。
(ql:quickload :lla)
外部ライブラリのパスを cl-user::*lla-configuration* に設定する。LLAをロードする以前に設定する必要があるので~/.sbclrc等の処理系の初期化ファイルに書いておくといい。ライブラリのパスは環境によって変わってるかもしれないのでlocate等で適当に調べる。
;; 普通のlibblasを使う場合 (defvar *lla-configuration* '(:libraries ("/usr/lib/libblas/libblas.so.3.0"))) ;; openblasを使う場合 (defvar *lla-configuration* '(:libraries ("/usr/lib/openblas-base/libblas.so.3"))) ;; intel MKLを使う場合 (defvar *lla-configuration* '(:libraries ("/usr/lib/x86_64-linux-gnu/libgomp.so.1" "/lib/x86_64-linux-gnu/libpthread.so.0" "/opt/intel/compilers_and_libraries_2016.1.150/linux/compiler/lib/intel64_lin/libiomp5.so" "/opt/intel/compilers_and_libraries_2016.1.150/linux/mkl/lib/intel64_lin/libmkl_rt.so" )))
行列積で各ライブラリの実行速度を比較
こちらの記事Python (NumPy) と Common Lisp (LLA) で行列積の実行速度を比較するのコードを使って実験する。
関数lla:mmで行列積が計算できる。行列はsingle-float型の2次元配列とする。ベクトルであっても2次元配列とする必要がある。
ランダムに初期化したN×N行列2つの掛け算を各Nについて100回実行した時間を計測する。環境はUbuntu14.04+Core i5 6440 (4コア)。
simpleは素のCommon Lispで型指定して最適化オプションなども付けて実装した行列積。simpleとblasはシングルスレッドなので遅い。メーカー謹製のMKLが一番速いと思っていたのだが、結果は何故かOpenBLASの方が速かった。
2016/2/19追記
こちらの記事Common Lispで高速行列演算に、上記の素のCommon Lispの行列積を、ループのインライン展開（アンローリング）やキャッシングによって更に最適化したもの(on-register-gemm)があったのでそれを図に追加してみた。
Lispbox for Windows (x86) バイナリ配布 2016年版
2014年版の内容もそろそろ古くなってきたかと思いアップデートした。Windows10 (64bit)で動作確認済み。
内容
構造体のcopy-関数は中身までコピーしてくれない
構造体のcopy-関数で中身のベクタもコピーされると思ったら大間違いという話。
(defstruct s (vec #(0000) :type simple-array)) (defparameter s1 (make-s)) (defparameter s2 (copy-s s1)) (eq s1 s2) ; => NIL (eq (s-vec s1) (s-vec s2)) ; => T
copy-sした上でベクタだけcopy-seqしたもので上書きすればいい。
(defun fullcopy-s (obj) (let ((clone (copy-s obj))) (setf (s-vec clone) (copy-seq (s-vec obj))) clone)) (defparameter s3 (fullcopy-s s1)) (eq s1 s3) ; => NIL (eq (s-vec s1) (s-vec s3)) ; => NIL
しかしいちいちこれを書くのはめんどさい。何かラッパーが欲しい
Common Lispで多層ニューラルネットを実装してみる
去年買ってはいたが積読状態だったこのディープラーニングの青い本。
- 作者: 岡谷貴之
- 出版社/メーカー: 講談社
- 発売日: 2015/04/08
- メディア: 単行本（ソフトカバー）
- この商品を含むブログ (13件) を見る
全ソースコードはGistに置いておく。(nn.lisp)
変数名などは本の中の数式と大体対応している。
データ構造
とりあえず一つのレイヤーの入力と出力の関係を記述するための構造体を用意する。ネットワーク全体はレイヤーのベクトルということにすればよろしい。
(defstruct layer in-dim ; 入力次元 out-dim ; 出力次元 (ユニット数) w-mat ; 出力次元 × 入力次元の重み配列 u-vec ; 総入力ベクトル z-vec ; 出力ベクトル delta-vec ; デルタ(誤差関数の総入力についての微分) activation-func ; 活性化関数 activation-func-diff) ; 活性化関数の導関数 (defstruct nn n-of-layers layer-vec learning-rate) ; 学習率
一つのニューロン(ユニット)に入ってくる入力に重みを掛けて和を取ったものが総入力uでありニューロンの興奮度合いを表す。uに活性化関数という普通はシグモイド型の関数をかませることで、ニューロンの出力はメリハリ(？)のあるものになる。
活性化関数はレイヤーごとに別々のものを使ってよいらしく、特に出力層では回帰なら恒等写像、二値分類ならロジスティック関数、多値分類ならソフトマックス関数と、問題に合わせて活性化関数を変えるらしい。
重みw-matをランダムに初期化してこれらの構造体のインスタンスを作る関数make-random-weight、make-random-layer、make-random-nnも定義しておく。
(defun make-random-weight (in-dim out-dim) (let ((w (make-array (list out-dim in-dim) :element-type 'double-float))) (loop for i from 0 to (1- out-dim) do (loop for j from 0 to (1- in-dim) do ;; initialize between -0.1 and 0.1 (setf (aref w i j) (- (random 0.2d0) 0.1d0)))) w)) (defun make-random-layer (in-dim out-dim activation-func activation-func-diff) (make-layer :in-dim in-dim :out-dim out-dim :w-mat (make-random-weight in-dim out-dim) :u-vec (make-array out-dim :element-type 'double-float :initial-element 0d0) :z-vec (make-array out-dim :element-type 'double-float :initial-element 0d0) :delta-vec (make-array out-dim :element-type 'double-float :initial-element 0d0) :activation-func activation-func :activation-func-diff activation-func-diff)) (defun make-random-nn (dimension-list activation-func-pair-list &optional (learning-rate 0.01d0)) (labels ((make-layers (product dimension-list activation-func-pair-list) (if (< (length dimension-list) 2) (nreverse product) (make-layers (cons (make-random-layer (car dimension-list) (cadr dimension-list) (caar activation-func-pair-list) (cadar activation-func-pair-list)) product) (cdr dimension-list) (cdr activation-func-pair-list))))) (make-nn :n-of-layers (1- (length dimension-list)) :layer-vec (apply #'vector (make-layers nil dimension-list activation-func-pair-list)) :learning-rate learning-rate)))
予測: 順伝播
まずは重みw-matは学習済みだと考えて、任意の入力から予測出力を計算するところをつくる。
総入力uは入力ベクトルとw-matの一つのユニットに対応する部分の線形結合を計算することによって、出力zはuに対して活性化関数をかませることで得られる。layer構造体のu-vecとz-vecを順に破壊的に更新することで一つのレイヤーの出力が計算できる。
(defun calc-u-vec (in-vec layer) (loop for j from 0 to (1- (layer-out-dim layer)) do (setf (aref (layer-u-vec layer) j) (loop for i from 0 to (1- (layer-in-dim layer)) summing (* (aref (layer-w-mat layer) j i) (aref in-vec i))))) (layer-u-vec layer)) (defun calc-z-vec (layer) (loop for i from 0 to (1- (layer-out-dim layer)) do (setf (aref (layer-z-vec layer) i) (funcall (layer-activation-func layer) (aref (layer-u-vec layer) i)))) (layer-z-vec layer))
あとは一つ前のレイヤーの出力を次のレイヤーの入力にして出力を計算するということを繰り返せばネットワーク全体の出力が計算できる。計算は入力層から出力層に向かって順方向に進むので順伝播と呼ぶ。
(defun forward (in-vec nn) (loop for i from 0 to (1- (nn-n-of-layers nn)) do (if (zerop i) (progn (calc-u-vec in-vec (aref (nn-layer-vec nn) i)) (calc-z-vec (aref (nn-layer-vec nn) i))) (progn (calc-u-vec (layer-z-vec (aref (nn-layer-vec nn) (1- i))) (aref (nn-layer-vec nn) i)) (calc-z-vec (aref (nn-layer-vec nn) i))))))
学習: 誤差逆伝播法
学習は、大枠としては普通の勾配法と一緒で、入力に対して順伝播で計算した予測と観測された教師信号の二乗誤差を重みで微分して勾配を得る。それから重みを勾配方向に少し更新する(この度合いを決めるのが学習率)。これを各データに対して繰り返すだけである。なので基本的にはオンライン処理になる。
この勾配を計算するために、誤差関数を総入力で微分した量(この本ではデルタと呼んでいる)が必要になるのだが、これが一つ後のレイヤーのデルタから計算できるというのがミソであり、出力層から入力層に向かって各ユニットのデルタを計算していく逆伝播のアルゴリズムが構成できる。これもlayer内のdelta-vecを破壊的に更新する。
(defun backward (train-vec nn) ;; calculate last layer's delta (let ((last-layer (aref (nn-layer-vec nn) (1- (nn-n-of-layers nn))))) (loop for j from 0 to (1- (layer-out-dim last-layer)) do (setf (aref (layer-delta-vec last-layer) j) (- (aref (layer-z-vec last-layer) j) (aref train-vec j))))) ;; calculate other deltas (loop for l from (- (nn-n-of-layers nn) 2) downto 0 do (let ((layer (aref (nn-layer-vec nn) l)) (next-layer (aref (nn-layer-vec nn) (1+ l)))) (loop for j from 0 to (1- (layer-in-dim next-layer)) do (setf (aref (layer-delta-vec layer) j) (* (funcall (layer-activation-func-diff layer) (aref (layer-u-vec layer) j)) (loop for k from 0 to (1- (layer-out-dim next-layer)) summing (* (aref (layer-delta-vec next-layer) k) (aref (layer-w-mat next-layer) k j)))))))))
最後に、forwardで計算したユニットの出力z-vec、backwardで計算したデルタdelta-vecを使って勾配を計算し、学習率をかけて重みw-matを更新する部分を書く。
(defun update (in-vec train-vec nn) (forward in-vec nn) (backward train-vec nn) ;; update first layer (let ((first-layer (aref (nn-layer-vec nn) 0))) (loop for i from 0 to (1- (layer-in-dim first-layer)) do (loop for j from 0 to (1- (layer-out-dim first-layer)) do (setf (aref (layer-w-mat first-layer) j i) (- (aref (layer-w-mat first-layer) j i) (* (nn-learning-rate nn) (aref in-vec i) (aref (layer-delta-vec first-layer) j))))))) ;; update other layer (loop for l from 1 to (1- (nn-n-of-layers nn)) do (let ((layer (aref (nn-layer-vec nn) l)) (prev-layer (aref (nn-layer-vec nn) (1- l)))) (loop for i from 0 to (1- (layer-in-dim layer)) do (loop for j from 0 to (1- (layer-out-dim layer)) do (setf (aref (layer-w-mat layer) j i) (- (aref (layer-w-mat layer) j i) (* (nn-learning-rate nn) (aref (layer-z-vec prev-layer) i) (aref (layer-delta-vec layer) j)))))))))
使い方
データを用意する
例によってサインカーブを近似する。入力ベクトルにはバイアスとして常に値が1になる次元を付けておく必要があることに注意(これをいつも忘れる)。
(defparameter input-data (loop repeat 100 collect (vector (- (random (* 2 pi)) pi) 1d0))) (defparameter train-data (mapcar (lambda (x) (vector (sin (aref x 0)))) input-data))
ネットワークを定義する
make-random-nnで重みw-matを初期化したネットワークを定義する。make-random-nnの第一引数は各層の次元数のリストであり、このリストの最初の要素が入力次元で、最後の要素が出力次元となる。第二引数は活性化関数と活性化関数の導関数の連想リストで、各レイヤーでどの活性化関数を使うかを指定する。第三引数は学習率であり、重みの更新幅に影響する。
とりあえず2層と3層の場合のネットワークの例を1つずつ示す。
;; 2層ネットワークの例(入力層を加えたら3層) (defparameter *nn* (make-random-nn '(2 50 1) ; 入力層2次元、隠れ層50次元、出力層1次元 (list (list #'RLF #'RLF-diff) ; 隠れ層の活性化関数: 正規化線形関数 (list #'identity #'one)) ; 出力層の活性化関数: 回帰問題なので恒等写像 0.05d0)) ; 学習率 ;; 3層ネットワークの例 (defparameter *nn* (make-random-nn '(2 20 10 1) (list (list #'logistic #'logistic-diff) ; 第1隠れ層の活性化関数をロジスティック関数にしてみる (list #'RLF #'RLF-diff) (list #'identity #'one)) 0.05d0))
1000サイクル学習する
(dotimes (i 1000) (mapc (lambda (in out) (update in out *nn*)) input-data train-data))
Common Lispで書く2個体分散遺伝的アルゴリズム
この記事はLisp Advent Calendar 2015 の25日目として書かれました。
遺伝的アルゴリズム(以下GA)は汎用的に使えるので、問題についての前提知識が何もないときにはとりあえずGAでやってみるということがよくある。
しかしCommon Lispのライブラリとしては以前紹介したGECOくらいしかない(と思う)。とはいえ実装するのはそんなに難しくないので、一から作ってみるのもいいかと思い、GA周りの比較的最近の話題を調べていると、とても実装が簡単で扱いやすそうな2個体分散遺伝的アルゴリズム(Dual DGA)(pdf)というモデルを見つけた。これを実装してみる。
分散遺伝的アルゴリズム(島モデル)
GAは計算が重い（特に評価部分）ので並列化して高速化しましょうということで、GAを分散化するモデルが色々とある。代表的なのは島モデルで、個体群を島と呼ばれるグループに分けて、島ごとに独立してGAをやるというものだ。各島は独立した計算になるので簡単に並列化できるし、高速化の度合いも大きい。
とはいえそのままだと島の中ですぐに多様性が失なわれて進化が収束してしまう。そこで、一定世代を経ると一定数の個体を他の島へ移住させることによって停滞してしまいがちな島社会に新風を吹き込んでやるということをする。これによって全体の多様性が保たれるので、初期収束を回避でき、普通にGAをやるよりもより良い解を見つけることが期待できる。速度と精度の両面で改善が見込めるうまい方法というわけだ。
一方で、島のサイズと移住間隔はメタパラメータとして与えてやらなければならない。ただでさえGAには交叉確率や突然変異確率、成績によって交配する個体を選択する方法などメタパラメータが多いので、これ以上メタパラメータが増えるのはうれしくない。
2個体分散遺伝的アルゴリズム
2個体分散遺伝的アルゴリズム(以下Dual DGA)は島モデルで島のサイズを2に固定した場合である。つまり総個体数/2の数だけ島がある。こうするとアルゴリズムが非常に単純化される。
島に2個体しかいないので、交配する相手は1人だけであり、選択方法を考える必要がない。交叉は必ず起き、遺伝子のどこか一箇所で必ず突然変異が起こるということにできる。問題を定義する部分以外には総個体数と移住間隔だけを決めればいいことになる。
論文からアルゴリズムの部分を抜粋すると、
1. 一つ前の移住から一定世代経過していた場合に，移住を行う．まず，2 つの個体のうち，ランダムに一方を選択し，そのコピーを他の島に送る．そして適合度の低い方の個体は，他の島から送られてきた個体に置き換えられる．
2. 2 つの個体を交叉させ，新しい 2 つの子個体を生成する．本論文では，一点交叉を用いている．この段階では，親個体も存在している．
3. 突然変異を行う．交叉で生成された 2 つの子個体をそれぞれ 1 ビット反転させるが，反転する点は，2 つの個体で 1 ビットずらす．これは，島内の 2 つの個体が同一になるのを防ぐためである．
4. 適合度の評価を行う．
5. 2 つの親個体と，2 つの子個体から，それぞれ適合度の高い方の個体を選び，次世代の 2 個体とする．
実装
レポジトリ https://github.com/masatoi/cl-dudga
上のアルゴリズムを200行くらいで書いてみたもの(Gist)。並列化のためにlparallelを使っているので事前にquickloadしておく必要がある。
lparallelのカーネルの設定は実行環境によって変える。
(defparameter *kernel* (make-kernel 4))
使い方
よくGAで使われるテスト関数Rastrigin関数を最小化する問題を解かせてみる。
評価関数の定義
まずは最大化（最小化）させたい評価関数を定義する。評価関数は一つの個体の遺伝子であるビットベクタを引数に取って、評価値を返す。
評価関数rastrign-evalを定義してみる。
(defun rastrign (x-list) (let ((n (length x-list))) (+ (* 10 n) (loop for xi in x-list summing (- (* xi xi) (* 10 (cos (* 2 pi xi)))))))) (defun bit-subseq->number (bit-vec start-index end-index) (let ((n (- end-index start-index))) (labels ((bin2dec (acc i n) (if (> i end-index) acc (if (= (bit bit-vec i) 1) (bin2dec (+ acc (expt 2 n)) (1+ i) (1- n)) (bin2dec acc (1+ i) (1- n)))))) (bin2dec 0 start-index n)))) (defun rastrign-eval (bit-vec) (let* ((dimension 20) (x-bit-length 10) (x-list (loop for i from 0 to (1- dimension) collect (bit-subseq->number bit-vec (* i x-bit-length) (1- (* (1+ i) x-bit-length)))))) (rastrign (mapcar (lambda (x) (* (- x 511.5) 0.01)) x-list))))
Rastrign関数は入力次元が2次元だとこういう形の関数になる。
実際には入力次元は20次元で実験する。各変数の定義域は[-5.12, 5.12]で、0.01刻みで10bitの二進数で表現する。なので20*10で200bitの二進数が遺伝子になる。
問題を定義する
次に、問題を定義する。problemとislandとindividualという構造体があり、problemはislandのベクタを持っており、islandはindividualのベクタを持っている。関数init-problemでランダムな個体で初期化されたproblemがつくられる。この問題では評価関数を最小化したいのでキーワード変数direction-of-optimizationにはminimizeを指定している（最大化の場合はmaximize)。
(defparameter prob (init-problem 300 ; 総個体数 200 ; 遺伝子長 5 ; 移住間隔 #'rastrign-eval ; 評価関数 (lambda (problem) ; 終了条件 (5000世代までいったら終了) (= (island-generation (aref (problem-islands problem) 0)) 5000)) :direction-of-optimization 'minimize))
GAを走らせる
最後に、probに対してrun-problemを実行することで計算が始まり、problem内の各個体が破壊的に変更されていく。
(time (run-problem prob)) ;; Evaluation took: ;; 3.393 seconds of real time ;; 12.470098 seconds of total run time (12.444631 user, 0.025467 system) ;; [ Run times consist of 0.064 seconds GC time, and 12.407 seconds non-GC time. ] ;; 367.52% CPU ;; 11,508,855,474 processor cycles ;; 3,461,615,296 bytes consed
1世代毎に全個体の平均評価値と最良の個体の評価値をプロットした図が以下になる。
ちゃんと最適化されていることが分かる。
TODO
遺伝子をグレイコード、実数ベクトルで表現してみる。交叉法を色々変える（今回実装したものは一点交叉）
ポール・グレアム Lispの起源 (The Roots of Lisp)
この記事はLisp Advent Calendar 2015 の10日目として書かれました。
- 作者: ポールグレアム,Paul Graham,川合史朗
- 出版社/メーカー: オーム社
- 発売日: 2005/01/01
- メディア: 単行本
- 購入: 109人 クリック: 4,884回
- この商品を含むブログ (582件) を見る
ハッカーと画家で有名なポール・グレアム氏が、ジョン・マッカーシーの1960年のLispのオリジナル論文をCommon Lispで解説した記事を書いていることに最近になって気付いた。
- The Roots of Lisp (www.paulgraham.com) (ソースコード)
かなり古い記事で(2002年)、内容はWikipediaの純LISPの項目と重なるところも多い。しかしながら現在でも非常に重要な記事なので、自分の勉強のためにも訳すことにした。別訳としてはlionfan氏による和訳があるが、こちらは前書き部分のみである。
この記事ではLispを一から定義して、evalを実装するところまでを解説している。ひとたびevalが実装されると、プログラムをあらゆる段階で評価できるようになるので、プログラムを変換してから評価したりできる。これは構文を新たに定義することに他ならないので、その上にどのような言語でも作ることができる。つまり、最も基本的な要素からメタプログラミングの土台を作るところまでを解説することが目的となっている。
この記事を理解できれば、他のどのような言語、計算機の上にも、自分自身でLispを実装できるようになる。一からプログラミング言語を作るという経験は非常に重要なので、ぜひ一度自分自身の手でLispを作ってみてほしい。
Lispの起源 (The Roots of Lisp)
1960年、ジョン・マッカーシーは驚くべき論文を発表した。彼はその論文で、ユークリッドが幾何学でやったのと同じようなことをプログラミングの世界でやった*1。 彼は少数の単純なオペレータと、関数のための記法を与え、そこからどのようにしてプログラミング言語全体を構築できるかを示した。彼はこの言語をリスト処理(LISt Processing)の意味を込めてLispと呼んだ。なぜなら彼のアイディアの中核の一つは、リストと呼ばれる単純なデータ構造をコードとデータの両方のために使うことだったからだ。
マッカーシーが何を発見したのかを理解することは、ただ単にコンピュータ史上のマイルストーンとして重要なだけでなく、今後プログラミングというものがどのように発展していくのかを考える上でのモデルとして価値がある。これまでのところ、本当の意味で混じりけなく一貫して継続しているプログラミングモデルは2系統あるように見える。すなわち、CモデルとLispモデルだ。これらはその間に沼地をはさんでそびえる2つの山のように見える。コンピュータがよりパワフルに成長するにつれて、発展を続ける新しい言語達は徐々にLispの山の方を目指して沼地の中を移動してきている。
過去20年間、新しくプログラミング言語をつくるときの人気のレシピはというと、基本的にはCモデルを採用して、Lispモデルから断片的に動的型付けやガベージコレクションのようなパーツを取ってきてはくっつけるというものだった。
この記事で私は、可能な限り単純な言葉でマッカーシーが何を発見したのかを説明するつもりでいる。ここで大事なのは、単に誰かが40年前に発見した興味深い理論的な結果について学ぶことではなく、今あるプログラミング言語達が将来的にどこを目指しているのかを示すことだ。Lispが奇妙に見えるのは(実際にはそれがLispの特質を決定づけているのだが)、「LispがLisp自身によって書ける」ということだ。これによってマッカーシーが何を示そうとしたのかを理解するために、我々は彼の数学的記法を実行可能なCommon Lispコードに置き換えながら、彼の足跡をたどっていくことにしよう。
7つの原始オペレータ
まず始めに式を定義する。式はアトムかリストのどちらかである。アトムは文字の並びである(例: foo)。リストは0個以上の式から成り、空白で区切られ、括弧でくくられている。ここでいくつか式の例を示す。
foo () (foo) (foo bar) (a b (c) d)
最後の式(a b (c) d)は4つの要素を持つリストだが、3番目の要素(c)はそれ自体が1つの要素を持つリストになっている。つまり、リストは入れ子にできる。
算術式1 + 1が値2を持つように、有効なLisp式もまた値を持つ。式eが値vを発生させるとき、それを「eがvを返す」と呼ぶ。式から値を取り出すことを評価と呼ぶ。例えば「式eを評価した結果、値vを得る」のように言う。
次のステップは、どのような種類の式がありうるか、そしてそれらがどのような値を返すかを定義することである。
ある式がリストであるとき、そのリストの最初の要素をオペレータと呼び、それ以降の要素を引数と呼ぶ。今から我々は7つの原始オペレータを定義していく。すなわち、quote、atom、eq、car、cdr、cons、そしてcondの7つだ。これらは数学でいう公理のような役割を持つ。
1. quote
(quote x)はxを返す。読みやすさのため、(quote x)を'xと略記することにする。
(quote a) ; => a 'a ; => a (quote (a b c)) ; => (a b c)
2. atom
(atom x)はxがアトムまたは空リストのときにアトムtを返し、それ以外では空リスト()を返す。Lispでは慣習的に、真理値を表す際は、アトムtを真とし、空リストを偽とする。
(atom 'a) ; => t (atom '(a b c)) ; => () (atom '()) ; => t
さて、今や我々はその引数が評価されるようなオペレータatomを手に入れたので、quoteが何のために使われるのかを示すことができる。その目的とは、リストをquoteすることでそのリストを評価から守ることである。atomのようなオペレータに、quoteされてないリストが引数として与えられた場合、それはコードとして扱われる。
(atom (atom 'a)) ; => t
逆に、quoteされたリストは評価から守られるため、単にリストとして扱われる。したがって、次の場合は2要素のリストになり、atomによるテストの結果は偽になる。
(atom '(atom 'a)) ; => ()
これはちょうど英語における引用符(quote)の使い方と似たようなものだ。単にCambridgeと書くとそれは「マサチューセッツ州の人口9万人の街」だが、引用符をつけて"Cambridge"と書くことで「9文字の英単語」になるのだ。つまり具体的な意味を持つ記述内容から表記上の文字列へと視点がシフトする。
他のプログラミング言語はquoteのような仕組みをほとんど持たないので、ちょっと異質な概念に思えるかもしれないが、それはLispが持つ最も独特な特徴の一つと密接に結びついている。その特徴とはすなわち、コードとデータがリストという同じデータ構造から出来ており、quoteオペレータがコードとデータを区別する方法であるということだ。
3. eq
比較のためのオペレータ。(eq x y)はxとyの値が同じアトムか、両方とも空リストであるときにtを返し、それ以外のときは()を返す。
(eq 'a 'a) ; => t (eq 'a 'b) ; => () (eq '() '()) ; => t
4. car
(car x)はxの値がリストであるとして、そのリストの最初の要素を返す。
(car '(a b c)) ; => a
5. cdr
(cdr x)はxの値がリストであるとして、そのリストの最初の要素を除いた残りの全ての要素を返す。
(cdr '(a b c)) ; => (b c)
6. cons
(cons x y)はyの値がリストであるとして、xの値にyの値の要素が続くようなリストを返す。
(cons 'a '(b c)) ; => (a b c) (cons 'a (cons 'b (cons 'c '()))) ; => (a b c) (car (cons 'a '(b c))) ; => a (cdr (cons 'a '(b c))) ; => (b c)
7. cond
条件分岐のためのオペレータ。(cond (p1 e1) ... (pn en))は以下のように評価される。p1からpnまでの式が順番に評価されていき、各条件式piの値がtになった時点で、それと対応するei式が評価され、その値がcond式全体の値として返される。
(cond ((eq 'a 'b) 'first) ((atom 'a) 'second)) ; => second
7つの原始オペレータのうち、quoteとcondを除く5つではその引数は常に評価される*2。 この種のオペレータを関数(function)と呼ぶ。
関数の表記法
次に、新たな関数を記述するための記法を定義する。関数は(lambda (p1 ... pn) e)のように記述される。ここでp1 ... pnはそれぞれパラメータと呼ばれるアトムであり、eは一つの式である。(各パラメータpiは式eの内部で参照できる)
((lambda (p1 ... pn) e) a1 ... an)
のように、関数をリストの最初の要素に置く式を関数呼び出しと呼ぶ。その値は以下のように計算される。まず各引数aiが評価される。次に、各パラメータpiの値が対応する引数aiの値に置き換えられ、それから関数本体部分eが評価される。
((lambda (x) (cons x '(b))) 'a) ; => (a b) ((lambda (x y) (cons x (cdr y))) 'z '(a b c)) ; => (z b c)
次の式のように、
(f a1 ... an)
式の第一要素が原始オペレータでないアトムfで、fの値が関数(lambda (p1 ... pn) e)であるなら、その式の値は
((lambda (p1 ... pn) e) a1 ... an)
の値となる。言い替えるなら、パラメータは式の中で引数として使われるのと同様にオペレータとしても使われうる。例えば、次の式におけるfは、'aをconsするという関数によって置き換えられ、オペレータとして'(b c)に対して作用する。
((lambda (f) (f '(b c))) (lambda (x) (cons 'a x))) ; => (a b c) ;; 訳注： このコードをCommon Lispで動かすためにはfがオペレータであることを宣言するためにfuncallが必要になる。 ;; Common Lispではパラメータとオペレータでは名前空間が異なるからだ。 ((lambda (f) (funcall f '(b c))) (lambda (x) (cons 'a x)))
これとはまた別に、関数が自分自身を参照できるようにする表記法があり、それによって再帰関数を定義するのが容易になる*3。 その記法は次のようなものだ。
(label f (lambda ( p1 ... pn ) e))
これは一つの関数を表しており、基本的には(lambda ( p1 ... pn ) e)のように振る舞うのだが、それに加えて、fが関数本体eの中であたかもこの関数のパラメータであるかのように現れたとき、そのfはこのlabel式自体へと評価されるという性質を持っている。
例えば、式x、アトムy、リストzを引数として取り、z中に(任意の深さで)出現するyをxに置き換えたリストを返す関数(subst x y z)を定義したいとする。
(subst 'm 'b '(a b (a b c) d)) ; => (a m (a m c) d)
この関数は以下のように表記できる。
(label subst (lambda (x y z) (cond ((atom z) (cond ((eq z y) x) ('t z))) ('t (cons (subst x y (car z)) (subst x y (cdr z)))))))
ここで、f = (label f (lambda ( p1 ... pn ) e))を次のように略記する。
(defun f ( p1 ... pn ) e)
そうすることで、substの定義はこうなる。(訳注： Common Lispにはこのlabel記法はないので、実際に動くのは以下のコードとなる。なお、substは組込み関数として既に定義されているため、subst.と置き換えた。ただし、以下で定義しているevalの内部ではlabel記法を扱える)
(defun subst. (x y z) (cond ((atom z) (cond ((eq z y) x) ('t z))) ('t (cons (subst. x y (car z)) (subst. x y (cdr z))))))
ちなみに、ここでcond式のデフォルト節をどうやって指定するかが出てきている。条件式が'tであるような節は常に成功する。従って、
(cond (x y) ('t z))
は他の言語でいうところのif x then y else zと等価である。
いくつかの補助関数
ここまでに我々は関数の表記法を手に入れたので、7つの原始オペレータを使って、新しい関数を定義していくことにする。まずは、いくつかの共通するパターンに省略形を導入することで便利になるだろう。そこで、carとcdrの組合せに対応する省略形としてcxxxrを使う。ここでxxxの部分はaかdの並びである。例えば、(cadr e)は(car (cdr e))の省略形であり、eの二番目の要素を返す。
(cadr '((a b) (c d) e)) ; => (c d) (caddr '((a b) (c d) e)) ; => e (cdar '((a b) (c d) e)) ; => (b)
また、(cons e1 ... (cons en '()) ...)の省略形として(list e1 ... en)を使う。
(cons 'a (cons 'b (cons 'c '()))) ; => (a b c) (list 'a 'b 'c) ; => (a b c)
さらにいくつか新しい関数を定義するが、原始オペレータと区別するため、またCommon Lisp組込みの関数との名前衝突を回避するために、これらの関数の名前の最後にはピリオド.を付けておくことにする。
1. null.
(null. x)はその引数xが空リストかどうかをテストする。
(defun null. (x) (eq x '())) (null. 'a) ; => () (null. '()) ; => t
2. and.
(and. x y)は両方の引数がtを返すときにtを返し、それ以外の場合は()を返す。
(defun and. (x y) (cond (x (cond (y 't) ('t '()))) ('t '()))) (and. (atom 'a) (eq 'a 'a)) ; => t (and. (atom 'a) (eq 'a 'b)) ; => ()
3. not.
(not. x)はその引数xが()を返すときにtを返し、その引数xがtを返すときに()を返す。
(defun not. (x) (cond (x '()) ('t 't))) (not (eq 'a 'a)) ; => () (not (eq 'a 'b)) ; => t
4. append.
(append. x y)は2つのリストを引数に取り、それらを連結したリストを返す。
(defun append. (x y) (cond ((null. x) y) ('t (cons (car x) (append. (cdr x) y))))) (append. '(a b) '(c d)) ; => (a b c d) (append. '() '(c d)) ; => (c d)
5. pair.
(pair. x y)は2つの同じ長さのリストを引数に取って、それぞれのリストの対応する位置にある要素のペアからなるリストを返す。
(defun pair. (x y) (cond ((and. (null. x) (null. y)) '()) ((and. (not. (atom x)) (not. (atom y))) (cons (list (car x) (car y)) (pair. (cdr x) (cdr y)))))) (pair. '(x y z) '(a b c)) ; => ((x a) (y b) (z c))
6. assoc.
(assoc. x y)は1つのアトムxと、pair.によって生成される形式のリストyを引数に取って、yの要素のリストのうち、第一要素がxであるようなリストの第二要素を値として返す。
(defun assoc. (x y) (cond ((eq (caar y) x) (cadar y)) ('t (assoc. x (cdr y))))) (assoc. 'x '((x a) (y b))) ; => a (assoc. 'x '((x new) (x a) (y b))) ; => new
サプライズ： evalの実装
ここまでのところで、リストを連結したり、一つの式を別の式に置き換えるなどの関数をエレガントな表記法で定義できるようになった。しかしそれが何になるというのだろうか？
ここで読者にサプライズがある。実はこの時点で既に我々の言語のためのインタプリタとして働く関数を書くことができる。インタプリタは任意のLisp式を引数に取ってその値を返す関数である。以下がその定義になる。
(defun eval. (e a) (cond ((atom e) (assoc. e a)) ((atom (car e)) (cond ((eq (car e) 'quote) (cadr e)) ((eq (car e) 'atom) (atom (eval. (cadr e) a))) ((eq (car e) 'eq) (eq (eval. (cadr e) a) (eval. (caddr e) a))) ((eq (car e) 'car) (car (eval. (cadr e) a))) ((eq (car e) 'cdr) (cdr (eval. (cadr e) a))) ((eq (car e) 'cons) (cons (eval. (cadr e) a) (eval. (caddr e) a))) ((eq (car e) 'cond) (evcon. (cdr e) a)) ('t (eval. (cons (assoc. (car e) a) (cdr e)) a)))) ((eq (caar e) 'label) (eval. (cons (caddar e) (cdr e)) (cons (list (cadar e) (car e)) a))) ((eq (caar e) 'lambda) (eval. (caddar e) (append. (pair. (cadar e) (evlis. (cdr e) a)) a))))) (defun evcon. (c a) (cond ((eval. (caar c) a) (eval. (cadar c) a)) ('t (evcon. (cdr c) a)))) (defun evlis. (m a) (cond ((null. m) '()) ('t (cons (eval. (car m) a) (evlis. (cdr m) a)))))
このeval.の定義は、これまでに見てきたどの定義よりも長いので、各部分がどのように働くかを詳細に見ていこう。
この関数は2つの引数を取る。eは評価される式で、aはこれまでに関数呼び出しにおけるパラメータとして出現してきたアトムとその値との対応のリストである。このリストは環境(environment)と呼ばれ、pair.によって生成される形式になっている。pair.やassoc.を書いたのは、これらのリストを構築し、探索する為だったのだ。
eval.の骨格は4つの節からなるcond式である。式をどのように評価するかは、その式がどんな種類かによる。このcond式の最初の節はアトムを取り扱う。もしeがアトムなら、環境の中からその値を見つけてくる。
(eval. 'x '((x a) (y b))) ; => a
2つ目の節は、aをアトムとして、(a ...)という形になっている式を取り扱うための別のcond式になっている。このcond式は全ての原始オペレータを含んでおり、それぞれに一つずつ対応する節がある。
(eval. '(eq 'a 'a) '()) ; => t (eval. '(cons x '(b c)) '((x a) (y b))) ; => (a b c)
これらのうちquoteを除く全てで、引数の値を見つけるために再帰的にeval.を呼んでいる。
最後の2つの節はより複雑になってくる。cond式を評価するために、補助関数evcon.を呼んでいる。evcon.はcond式の節のリストcを再帰的に調べていき、節の最初の要素(条件式)がtを返すような節を見つけ、二番目の要素の値を返す。
(eval. '(cond ((atom x) 'atom) ('t 'list)) '((x '(a b)))) ; => list
eval.の二番目の節の最後の部分(つまり2つ目のcond式のデフォルト節)は、パラメータとして渡されている関数の呼び出しを取り扱う。これは、(a ...) という形のアトムaをその値(lambda式かlabel式であるべき)で置き換えて、その結果を評価することによって動作する。従って、
(eval. '(f '(b c)) '((f (lambda (x) (cons 'a x)))))
は
(eval. '((lambda (x) (cons 'a x)) '(b c)) '((f (lambda (x) (cons 'a x)))))
となって、(a b c)を返す。
eval.の最後の二つの節は、式の最初の要素にlambda式またはlabel式を直接置くような関数呼び出しを取り扱う。label式の場合、まず関数名と関数本体のリストを環境に追加し、それからlabel式内部のlambda式によってこのlabel式自身を置き換え、その式に対してeval.を呼ぶことにより評価される。つまり、
(eval. '((label firstatom (lambda (x) (cond ((atom x) x) ('t (firstatom (car x)))))) y) '((y ((a b) (c d)))))
は
(eval. '((lambda (x) (cond ((atom x) x) ('t (firstatom (car x))))) y) '((firstatom (label firstatom (lambda (x) (cond ((atom x) x) ('t (firstatom (car x))))))) (y ((a b) (c d)))))
になり、実際にはaを返す。
最後に、((lambda (p1 ... pn) e) a1 ... an)の形の式は次のようにして評価される。まず引数(a1 ... an)の値のリスト(v1 ... vn)を得るためにevlis.を呼び、それから(p1 v1) ... (pn vn)を環境の先頭に追加してeを評価する。だから、
(eval. '((lambda (x y) (cons x (cdr y))) 'a '(b c d)) '())
は
(eval. '(cons x (cdr y)) '((x a) (y (b c d))))
になり、実際には(a c d)を返す。
まとめ
今や我々はevalがどのようにして動くかを理解したので、立ち戻ってそれが何を意味しているのかを考えてみよう。ここで我々が得たものは、非常にエレガントな計算のモデルである。我々は、ただquote、atom、eq、car、cdr、cons、condのみを使って、eval.を定義することができた。そしてそれは実際に我々の言語を実装するものであり、それを使って我々が欲するどのような追加の関数も定義することができる。
もちろん、既に様々な計算のモデルは存在している。最も有名なものはチューリングマシンだが、チューリングマシンのプログラムはあまり読みやすくなるようには考えられていない。もしあなたがアルゴリズムを記述するための言語を欲しているなら、記述をより抽象化するためのなにかしらの手段が欲しいと思うだろう。そしてそれがマッカーシーがLispを定義した狙いの一つなのだ。
1960年に定義された言語に足りないものは多い。まず副作用が無いし、副作用を使うときに便利な順次実行の仕組み(つまりCommon Lispのprogn、Schemeのbegin)もない。実用的な数値もないし*4、ダイナミックスコープもない。しかしこれらの制限は、驚くほどわずかなコードを追加するだけで取り除くことができる。SteeleとSussmanはそれをどうやるかを"The Art of the Interpreter"と呼ばれる有名な論文で示した*5。
もしあなたがマッカーシーのevalを理解したのなら、単にプログラミング言語史の一つのステージを理解したという以上のものを理解していることになる。これらのアイディアは今日のLispにおいても意味論的な中核であるからだ。マッカーシーの元論文を研究することは、ある意味で、Lispが本当は何であるのかを我々に示してくれる。Lispはマッカーシーがデザインしたものというよりは発見したものという方が近い。Lispは本来、人工知能やラピットプロトタイピング(あるいは他のそのようなレベルのタスク)のための言語というわけではない。Lispとは、計算を公理化しようと試みた結果として得られる何かなのである。
長い時間をかけて、平均的なプログラマによって使われる言語は継続してLispに近付いてきている。evalを理解することによって、将来の主流な計算のモデルがどうなっていくかを推測することができるだろう。
*1:"Recursive Functions of Symbolic Expressions and Their Computation by Machine, PartI." Communications of the ACM 3:4, April 1960, pp. 184-195.
*2:quoteとcondで始まる式は異なる評価のされ方をする。quote式が評価されると、その引数は評価されずにそのままquote式全体の値として返される。また、正しいcond式では、条件にマッチした部分式だけが評価され、cond式全体の値として返される。
*3:論理的には、このために新しい表記法を定義する必要はない。ここまでに定義した記法を使って、Yコンビネータと呼ばれる関数を扱う関数を使うことによって再帰関数を定義することができる。マッカーシーはこの論文を書いた時点でYコンビネータについては知らなかったかもしれないが、なんであれlabel記法はより読みやすい。
*4:マッカーシーの1960年版Lispでも、例えば数値nを表現するのにn個のアトムを持つリストを使うなどすれば可能ではある
*5: Guy Lewis Steele, Jr. and Gerald Jay Sussman, "The Art of the Interpreter, or the Modularity Complex - Parts Zero, One, and Two," MIT AI Lab Memo 453, May 1978.
書籍「オンライン機械学習」を買ったのでCommon Lispで実装してみた。(AROW編)
前回の記事ではパーセプトロンと線形SVMを実装したが、より高度な手法として、CW、AROW、SCWといったものがあるらしい。
パーセプトロン等では重みベクトルを直接更新していたのに対して、CWでは重みベクトルが正規分布に従って分布していると仮定して、その正規分布の平均や分散を更新していく。この平均と分散はそれぞれ重みベクトルと特徴量ごとの信頼度に対応し、これがCW(Confidence Weighted)の名前の由来となっている。更新にあたっては、誤分類の確率を一定値以下にするという制約を付けた上で、推定された正規分布とのカルバックライブラーダイバージェンス(分布間の距離的なもの)を最小化するような正規分布の平均と分散を求める。
今回実装するAROWやSCWはCWの派生で、CWがノイズが弱いという問題点を解決した手法である。
AROWやSCWについては解説した記事があるのでそちらを。
- AROWやSCWのまとめ http://kazoo04.hatenablog.com/entry/2012/12/20/000000
- SCWの論文のまとめ http://d.hatena.ne.jp/kisa12012/20120625/1340616659
- AROWの論文 http://webee.technion.ac.il/people/koby/publications/arow_nips09.pdf
- SCWの論文 http://icml.cc/2012/papers/86.pdf
精度はSCW > AROWだけど、SCWの論文中の比較表を見るにそれほど大した差ではないし、AROWの方が優れている場合もある。SCWはハイパーパラメータが2つでグリッドサーチしなければならないのに対して、AROWはハイパーパラメータが1つで扱いやすい。
AROWの疑似コード
実装上の注意点としては、共分散行列の更新は行列演算なので素直にやると計算量がO(m^3)になるが、対角成分以外を0として対角行列に近似することでO(m)にできるということがある。実際には、対角行列は行列ではなくベクトルとして実装し、対角行列とベクトルの乗算をする関数diagonal-matrix-multiplicationを定義しておく。
すると上記の疑似コードの対応する1ステップの更新部分はこうなる。
;; AROW、入力ベクトルに1次元足すバージョン(バイアスを分けない) ;; muとsigmaが破壊的に更新される。gammaがハイパーパラメータ (defun train-arow-1step (input mu sigma gamma training-label tmp-vec1 tmp-vec2) (let ((loss (- 1d0 (* training-label (inner-product mu input))))) ; ヒンジ損失が0より大きいときに更新 (if (> loss 0d0) (let* ((beta (/ 1d0 (+ (inner-product (diagonal-matrix-multiplication sigma input tmp-vec1) input) gamma))) (alpha (* loss beta))) ;; muの更新 ;; betaの計算の時点で、tmp-vec1には Sigma_{t-1} x_t の結果が入っている ;; muの更新差分をtmp-vec2に入れる (v-scale tmp-vec1 (* alpha training-label) tmp-vec2) (v+ mu tmp-vec2 mu) ;; sigmaの更新 ;; \Sigma_{t-1} x_t x_t^T \Sigma_{t-1} を対角行列に近似したものをtmp-vec1に入れる (diagonal-matrix-multiplication tmp-vec1 tmp-vec1 tmp-vec1) ;; betaをかけてtmp-vec1を更新 (v-scale tmp-vec1 beta tmp-vec1) ;; sigmaを更新 (v- sigma tmp-vec1 sigma))) (values mu sigma)))
このままだとバイアスに対応する平均と分散がmuとsigmaに含まれているため、データセットの各入力ベクトルに1次元、定数要素を追加する前処理が必要になりなんか嫌である。バイアスの平均と分散の更新は分けて、2要素のベクトルmu0-sigma0-vecに格納することにすると、上の関数はこうなる。
(defun train-arow-1step (input mu sigma mu0-sigma0-vec gamma training-label tmp-vec1 tmp-vec2) (let* ((mu0 (aref mu0-sigma0-vec 0)) (sigma0 (aref mu0-sigma0-vec 1)) (loss (- 1d0 (* training-label (f input mu mu0))))) ; ヒンジ損失が0より大きいときに更新 (if (> loss 0d0) (let* ((beta (/ 1d0 (+ sigma0 (inner-product (diagonal-matrix-multiplication sigma input tmp-vec1) input) gamma))) (alpha (* loss beta))) ;; muの更新 ;; betaの計算の時点で、tmp-vec1には Sigma_{t-1} x_t の結果が入っている ;; muの更新差分をtmp-vec2に入れる (v-scale tmp-vec1 (* alpha training-label) tmp-vec2) (v+ mu tmp-vec2 mu) ;; mu0の更新 (setf (aref mu0-sigma0-vec 0) (+ mu0 (* alpha sigma0 training-label))) ;; sigmaの更新 ;; Sigma_{t-1} x_t x_t^T Sigma_{t-1} を対角行列に近似したものをtmp-vec1に入れる (diagonal-matrix-multiplication tmp-vec1 tmp-vec1 tmp-vec1) ;; betaをかけてtmp-vec1を更新 (v-scale tmp-vec1 beta tmp-vec1) ;; sigmaを更新 (v- sigma tmp-vec1 sigma) ;; sigma0の更新 (setf (aref mu0-sigma0-vec 1) (- sigma0 (* beta sigma0 sigma0))))) (values mu sigma mu0-sigma0-vec)))
あとはデータセットにこれを適用する部分を書くだけである。
(defun train-arow-all (training-data mu sigma mu0-sigma0-vec gamma tmp-vec1 tmp-vec2) (loop for datum in training-data do (train-arow-1step (cdr datum) mu sigma mu0-sigma0-vec gamma (car datum) tmp-vec1 tmp-vec2)) (values mu sigma mu0-sigma0-vec)) (defun train-arow (training-data gamma) (let* ((dim (length (cdar training-data))) (mu (make-dvec dim 0d0)) (sigma (make-dvec dim 1d0)) (mu0-sigma0-vec (make-array 2 :element-type 'double-float :initial-contents '(0d0 1d0))) (tmp-vec1 (make-dvec dim 0d0)) (tmp-vec2 (make-dvec dim 0d0))) (train-arow-all training-data mu sigma mu0-sigma0-vec gamma tmp-vec1 tmp-vec2)))
前回同様にして読み込んだデータセットで試してみると
(time (multiple-value-bind (mu sigma mu0-sigma0-vec) (train-arow a1a-train 10d0) (declare (ignore sigma)) (test a1a-test mu (aref mu0-sigma0-vec 0)))) ; Accuracy: 84.461815%, Correct: 26146, Total: 30956 ; Evaluation took: ; 0.013 seconds of real time ; 0.013124 seconds of total run time (0.013124 user, 0.000000 system) ; 100.00% CPU ; 42,935,095 processor cycles ; 1,900,544 bytes consed (time (multiple-value-bind (mu sigma mu0-sigma0-vec) (train-arow a9a-train 10d0) (declare (ignore sigma)) (test a9a-test mu (aref mu0-sigma0-vec 0)))) ; Accuracy: 84.94564%, Correct: 13830, Total: 16281 ; Evaluation took: ; 0.032 seconds of real time ; 0.032499 seconds of total run time (0.032499 user, 0.000000 system) ; 100.00% CPU ; 107,560,608 processor cycles ; 9,076,736 bytes consed
多少計算時間が増えるものの、適当に設定したgammaでも良好な精度を出すことが分かる。線形SVMよりハイパーパラメータも減って精度もいいので常用できそうである。
今後は、SCW、平均化確率的勾配降下法(ASGD)、疎ベクトルの演算あたりをやるかもしれない。
