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

遺伝子をグレイコード、実数ベクトルで表現してみる。交叉法を色々変える(今回実装したものは一点交叉)