Pythonでの実数値遺伝的アルゴリズム(Real-coded genetic algorithm)の実装

実数値遺伝的アルゴリズム(Real-coded genetic algorithm)について勉強したのでその記録です。
以下、実数値遺伝的アルゴリズムをRCGAとしています。

以下の検証に関するコードはgithubにあげてあります。
github.com

ちなみにこちらのページで、今回実装したRCGAを強化学習に適用しています。
st1990.hatenablog.com

1. 本記事の概要

RCGAの概要
RCGA(JGG-AREX)のpython実装
RCGAのテスト

2. RCGAの概要

遺伝的アルゴリズム

最適化問題とは、とある関数f(x)について、f(x)が最小となるxを求める問題のことです。
例えば、ディープラーニングのパラメータの学習も最適化問題です。学習データと推定結果の差異に関する関数L(コスト関数と呼ばれています)を最小とするパラメータを求めています。

最適化問題を解くアルゴリズムとして代表的なものに勾配法があります。
勾配法は最適化したい関数f(x)の勾配(微分値)を使って、f(x)を最小とするようなxを求めます。
先述のディープラーニングの学習も勾配法が使われています。

遺伝的アルゴリズムでは、最適化問題の解を遺伝子に見立て、遺伝子の世代交代や組み換えを繰り返すことでより良い遺伝子（＝最適化問題の良い解）を求めます。
具体的には以下のような操作をします。

① 解をベクトルで表現する（このベクトルが遺伝子と呼ばれる）。ただし、解は整数で表現できるものに限ります。例)例えば、a,b,cの3要素で値が決まる関数f(a,b,c)であれば、遺伝子gはg=[a,b,c]で表現される。
② 遺伝子を複数用意する(現世代群)。
③ 現世代群から親となる遺伝子を選択する。
④ 選択した遺伝子を交叉（配合みたいなイメージ）させ、子の遺伝子を作成する。
⑤ 子の遺伝子の良さ（最小にしたい関数f(x)の値がどれだけ小さいか）を評価し、目標基準以上となれば終了。
⑥ 現世代群、親遺伝子、子遺伝子から次世代群となる遺伝子を選び、それを現世代群とする。
⑦ ③に戻る。
ただし、世代交代方法や交叉方法にはいろいろな種類があり、問題ごとに適した手法を見極めた方が良いとか。

遺伝的アルゴリズムのメリットとデメリットは以下のとおりです。

メリット：最適化したい関数f(x)の勾配が必要ない（=定式化や微分の計算が困難な問題でも取り扱える、関数f(x)をブラックボックスとして取り扱える）。
デメリット：勾配を使わないので収束が遅い。

Google先生に聞くと詳しくてわかりやすい説明がたくさん出てきます。

実数値遺伝的アルゴリズム(RCGA)

普通の遺伝的アルゴリズムは実数値を取り扱えませんでしたが、RCGAは実数値を取り扱えるようになっています。
RCGAも遺伝的アルゴリズムと大まかな流れは変わりませんが交叉の方法がだいぶ違います。
今回は世代交代方法としてJGG、交叉方法としてAREXを実装したので、これらの手法の概要イメージを載せておきます。

f:id:st1990:20190420143305p:plain — JGGの世代交代イメージ

f:id:st1990:20190420143353p:plain — RCGA(AREX)の交叉イメージ

詳しいことは以下のサイトを参照ください。

AREXの原論文(ハイパーパラメータの推奨値のってます) https://www.jstage.jst.go.jp/article/tjsai/24/6/24_6_446/_pdf
AREXを使っている論文 http://t2r2.star.titech.ac.jp/rrws/file/CTT100683247/ATD100000413/
実数値遺伝的アルゴリズム実装ざっくり説明 - Qiita

3. RCGAの実装

import numpy as np

class RealCodecGA_JGG_AREX:
    '''
    実数値遺伝的アルゴリズム（生存選択モデルJGG、交叉SPX）

    初期化
    交叉する個体の選択
    交叉
    残す個体の選択
    '''
    def __init__(self, gene_num, 
                 evaluation_func, 
                 initial_min, initial_max, 
                 population=None, 
                 crossover_num=None, child_num=None, 
                 initial_expantion_rate=None, learning_rate=None, 
                 seed=None):
        '''
        recommended value of parameter are following.
         crossover_num = gene_num + 1
         child_num = 4 * gene_num (~4 * crossover_num)
         initial_expantion_rate = 1.0
         learning_rate = 1/5/gene_num
        '''

        #乱数シード
        np.random.seed(seed=seed)

        # 遺伝子情報
        self.gene_num = gene_num

        # 個体の評価関数
        # evaluation function defined as following.
        #  func(genes):
        #    return evaluations
        self.evaluation_func = evaluation_func

        #遺伝子の初期値範囲
        self.initial_min = initial_min
        self.initial_max = initial_max

        # 全個体数
        # 交叉する個体数
        if (population is not None) and (crossover_num is not None):
            self.population = population
            self.crossed_num = crossover_num
        elif (population is None) and (crossover_num is not None):
            self.crossed_num = crossover_num
            self.population = self.crossed_num
        elif (population is not None) and (crossover_num is None):
            self.population = population
            self.crossed_num = np.minimum(gene_num+1, population)
        elif (population is None) and (crossover_num is None):
            self.crossed_num = gene_num+1
            self.population = self.crossed_num

        # 交叉して生成する子個体数
        self.child_num = child_num if child_num is not None else 4*self.crossed_num

        # 拡張率
        self.expantion_rate = initial_expantion_rate if initial_expantion_rate is not None else 1.0
        # 拡張率の学習率
        self.learning_rate = learning_rate if learning_rate is not None else 1.0/5.0 / gene_num

        #
        self.genes = None # 候補の遺伝子
        self.evals = None # 候補の遺伝子の評価値
        self.best_gene = None # 現時点で評価値が一番高い遺伝子
        self.best_evaluation = None # 現時点で一番高い評価値
        self.last_gene = None # 最終時点での遺伝子
        
        #
        self.__initialize_genes()

        return

    def __initialize_genes(self):
        # 遺伝子の初期値設定
        # genes = [[gene_0], [gene_1], ... ,[gene_population]]
        self.genes = self.initial_min + (self.initial_max - self.initial_min) * np.random.rand(self.population, self.gene_num)
        
        # 遺伝子の評価値
        self.evals = self.evaluation_func(self.genes)
        
        # min
        min_idx = np.argmin(self.evals)
        self.best_evaluation = self.evals[min_idx]
        self.best_gene = self.genes[min_idx].copy()

        return

    def generation_step(self):
        # 個体数
        pop = len(self.genes)

        # 交叉する個体の選択
        # 交叉する個体のインデックス
        crossed_idx = self.__random_select(pop, self.crossed_num)

        # 交叉
        child_genes, rand_mtrx = self.__arex_crossover(self.genes[crossed_idx], self.evals[crossed_idx], self.child_num, self.expantion_rate)

        # 残す個体の選択
        survive_genes, survive_evals, survive_idx = self.__ranking_survival(child_genes, self.evaluation_func, survive_num=self.crossed_num)
        # 親個体と変更
        self.genes[crossed_idx] = survive_genes.copy()
        self.evals[crossed_idx] = survive_evals

        # 拡張率の更新
        self.expantion_rate = self.__update_arex_expantion_rate(self.expantion_rate, self.learning_rate, rand_mtrx[survive_idx], crss_num=self.crossed_num)
        #print(self.expantion_rate)

        min_idx = np.argmin(survive_evals)
        # 最終遺伝子の更新
        self.last_gene = survive_genes[min_idx].copy()
        # 生き残ったうちの最高評価
        best_survive_evals = survive_evals[min_idx]
        # 最高評価の更新
        if self.best_evaluation > survive_evals[min_idx]:
            self.best_evaluation = survive_evals[min_idx]
            self.best_gene = survive_genes[min_idx].copy()
        
        return best_survive_evals

    @staticmethod
    def __arex_crossover(genes, evals, child_num, expantion_rate):
        # 交叉される遺伝子数
        crss_num = len(genes)

        # 評価が高い順に並べたインデックス
        sorted_idx = np.argsort(evals)

        # 荷重重心
        w = 2.0 * (crss_num + 1.0 - np.arange(1, crss_num+1)) / (crss_num * (crss_num + 1.0))
        wG = np.dot(w[np.newaxis,:], genes[sorted_idx])

        # 重心
        G = np.average(genes, axis=0)

        # 乱数
        rnd_mtrx = np.random.normal(loc=0.0, scale=np.sqrt(1/(crss_num-1)), size=(child_num, crss_num))
        
        # 子個体
        child_genes = wG + expantion_rate * np.dot(rnd_mtrx, genes - G)

        return child_genes, rnd_mtrx

    @staticmethod
    def __update_arex_expantion_rate(expantion_rate, learning_rate, survive_rand_mtrx, crss_num=None):
        survive_num = len(survive_rand_mtrx)
        #
        crss_num_ = crss_num if crss_num is not None else survive_num
        #
        ave_r = np.average(survive_rand_mtrx, axis=0)
        L_cdp = expantion_rate**2 * (survive_num - 1) * (np.sum(np.square(ave_r)) - (np.sum(ave_r))**2 / survive_num)
        L_ave = expantion_rate**2 * (1/(crss_num_ - 1)) * (survive_num - 1)**2 / survive_num
        #
        new_expantion_rate = np.maximum(1.0, expantion_rate * np.sqrt((1.0 - learning_rate) + learning_rate * L_cdp / L_ave))

        return new_expantion_rate

    @staticmethod
    def __ranking_survival(genes, evaluation_func, survive_num):
        evals = evaluation_func(genes)
        survive_idx = np.argpartition(evals, survive_num)[:survive_num]
        #
        survive_genes = genes[survive_idx].copy()
        survive_eval = evals[survive_idx]
        return survive_genes, survive_eval, survive_idx

class Sample:
    @staticmethod
    def sample_RealCodecGA_JGG_AREX():
        # define evaluation func
        def minus_l2(ndar):
            ev = - np.sum(np.square(ndar), axis=1)
            return ev
        # rcga class
        rcga = RealCodecGA_JGG_AREX(gene_num=2, 
                                    evaluation_func=minus_l2, 
                                    initial_min=-1, initial_max=1, 
                                    population=10, 
                                    crossover_num=3, child_num=6, 
                                    initial_expantion_rate=1.0, learning_rate=0.05, 
                                    seed=1000)
        # optimize
        for i in range(100):
            rcga.generation_step()
            print('{0}'.format(i+1))
            print(' best evals : {0}'.format(rcga.best_evaluation))
            print(' best genes : {0}'.format(rcga.best_gene))
        return

4. RCGAのテスト

RCGAのテストをしました。
テストケースは①二次関数の最小化、②重回帰係数の計算の2つです。

テスト①：二次関数の最小化

二次関数を最小化する問題をRCGAで解きます。

最適化関数は以下のとおり。
$f(a _ {1}, a _ {2}, ..., a _ {100}) = (a _ {1}^{2}+a _ {2}^{2}+...+a _ {100}^{2})/100$
解はすべてのa=0です。

RCGAで解くため、遺伝子を[a _ {1}, a _ {2}, ..., a _ {100}]で表しました。
aの初期値は[-1, 1]の範囲の一様乱数から生成しました。
ハイパーパラメータは世代の個体数を400として、あとは原論文の推奨値としました。
（親個体数101、子個体数404=4*101）

収束の様子を下図に示します。
a=0が解なので、f=0となれば解が求まったと言えます。

f:id:st1990:20190420151009p:plain — 二次関数最小化の収束の様子(n=100、縦軸は現世代までの最小のfの値、横軸は世代数)

…おそ！！覚悟はしていましたが、勾配法と比べると遅いです。この計算に2秒はかかってます。
これだけ計算して、評価値が10^-4か…xがやっと0.01ぐらいになりました。。
各世代404回fを計算しているので、世代数500までに合計202000回fを計算したことになります。fの計算が重いようなケースだとやばいです。

でもグリッドサーチよりはだいぶいいですね。[-1,1]の範囲を10分割したとして変数100個なので、計算終わるまでに10¹⁰⁰回fを計算しないといけないです。

テスト②：重回帰係数の計算

重回帰問題の係数をRCGAで求めます。

問題設定は以下のとおりです。

目的変数1個、説明変数100個、サンプル数5000個。
目的変数を $y _ {i}$ 、説明変数を $x _ {i,1},x _ {i,2},...,x _ {i,100}$ としました（i=1～5000）。
$y _ {i}$ は $y _ {i}=b^{ref}+a _ {1}^{ref}x _ {i,1}+a _ {2}^{ref}x _ {i,2}+...+a _ {100}^{ref}x _ {i,100}+ε _ {i}$ で計算したものを使いました。
$b^{ref}, a _ {1}^{ref}, a _ {2}^{ref}, ..., a _ {100}^{ref}$ は[-1, 1]の一様乱数からサンプリングしました。
$ε _ {i}$ はノイズで、平均0、標準偏差0.01の正規乱数からサンプリングしました。

回帰の係数を $b, a _ {1}, a _ {2}, ..., a _ {100}$ とすると、最適化関数を次式で表せます。
$f(b, a _ {1}, a _ {2}, ..., a _ {100}) = \sum_{i=0}^{100} (y _ {i} - y _ {i}^{pre})^{2} / 100$
$y _ {i}^{pre}$ は推定値で、 $y _ {i}^{pre}=b+a _ {1}x _ {i,1}+a _ {2}x _ {i,2}+...+a _ {100}x _ {i,100}$ です。

RCGAで解くため、遺伝子を[b, a _ {1}, a _ {2}, ..., a _ {100}]で表しました。
b, aの初期値は[-1, 1]の範囲の一様乱数から生成しました。
ハイパーパラメータは世代の個体数を1000として、あとは原論文の推奨値としました。
（親個体数101、子個体数404=4*101）

収束の様子を下図に示します。
$y _ {i}$ には標準偏差0.01のノイズがのっているので、概ねf=0.0001となれば解が求まったと言えます。

f:id:st1990:20190421210007p:plain — 重回帰係数の計算(n=100, σ=0.01、縦軸は現世代までの最小のfの値、横軸は世代数)

これもかなり遅いですが、なんとか解を求められました。この計算に10秒以上かかっています。
勾配法であれば一瞬です。
GAは最適化関数をブラックボックスで取り扱えるという強力なメリットがある反面、計算速度は遅いですね。勾配法を使える時はおとなしく勾配法を使った方がいいです。

遺伝子の要素数がかなり多くなる時は計算時間が膨大になるので、世代個体数や親個体数などのハイパーパラメータをちゃんと調整してやらんと厳しそうです。

5. まとめ

実数値遺伝的アルゴリズム(RCGA)の概要を説明しました。
RCGAをpythonで実装しました。
RCGAを二次関数の最小化、重回帰係数の計算に適用し、性能をテストしました。
- どちらの問題でも解を求めることができました。
- 勾配法と比べると計算速度は圧倒的に遅いです。勾配法が使えるときは勾配法を使いましょう。使えないときはRCGAは良い選択肢になると思います。
- 遺伝子の要素数が多くなるような問題では、ハイパーパラメータをしっかり調整してやらないと計算時間がめちゃくちゃかかりそうです。

statsuのblog

愛知のデータサイエンティスト。自分の活動記録。主に機械学習やその周辺に技術について学んだことを記録していく予定。