分布推定アルゴリズムを学ぶ (1) ～概要～
分布推定アルゴリズムを学ぶ (2) ～PBIL～
分布推定アルゴリズムを学ぶ (3) ～UMDA～
分布推定アルゴリズムを学ぶ (4) ～CGA～
分布推定アルゴリズムを学ぶ (5) ～ECGA～の続き

Bayesian Optimization Algorithm(BOA)とは

目的関数の変数をベイジアンネットワークを用いて表現することで，変数間に依存関係を持つ問題に対処した分布推定アルゴリズム手法．
ECGAと同様に，離散空間全体を対象とした最適化手法．

概要

BOAでは，確率モデルとしてベイジアンネットワーク(Bayesian Network; BN)を用い，変数間の依存関係を表現しながら最適化を行う．
ベイジアンネットワークは，ネットワーク構造と各ノードの条件付き確率から構成され，BOAではこれらを集団から推定する必要がある．

以下では，ベイジアンネットワークの概要，ベイジアンネットワークのネットワーク構造の構築，各ノードの条件付き確率の推定，ベイジアンネットワークからの個体の生成について説明し，BOAのアルゴリズムについてまとめる．

ベイジアンネットワーク

ベイジアンネットワークは多項データをモデリングする際によく用いられるモデルである．
ベイジアンネットワークは，変数間の依存関係を木構造によって表現する．
以降，$i$次元目の変数($i$番目のノード)を$X_i$と表記する．

ベイジアンネットワークは，非巡回有向グラフを仮定しており，以下の式によって表される．

P(X)=\Pi_{i=1}^{D}P(X_i|\Pi_{X_i}),

ここで，$X=(X_1,...,X_D)$は$D$次元ベクトル，$\Pi_{X_i}$はノード$X_i$の親ノードの集合である．

以下に，3次元ビット列の場合のベイジアンネットワークの例を示す．
親ノードの集合はそれぞれ$\Pi_{X_1}=\{\}$，$\Pi_{X_2}=\{X_1\}$，$\Pi_{X_3}=\{X_1,X_2\}$であり，各ノードは親ノードに対する条件付き確率を持つ．

ネットワーク構造の構築

ベイジアンネットワークのネットワーク構造$B$はK2アルゴリズムを用いて構築する．
K2アルゴリズムでは，任意のノード$X_i$とその親ノードの集合$\Pi_{X_i}$の間の依存度の高さを任意の評価指標Score($X_i,\Pi_{X_i}$)を用いて計算し，スコアを最大化することで最も集団にフィットしたネットワーク構造を探索する．
K2アルゴリズムの流れは以下の通りである．
ここで示すアルゴリズムは，オリジナルのK2アルゴリズムとは若干異なる点に注意されたい．

各ノードの親ノードの集合を空集合で初期化．つまり，$\Pi_{X_i}\leftarrow\emptyset$．
任意のノード$X_i$と$X_j$が依存関係を持っていると仮定，すなわち$\Pi_{X_i}^j\leftarrow\Pi_{X_i}\cup\{X_j\}$として，全ての組に対してScore($X_i,\Pi_{X_i}^j$)を計算．
全ての組の中で最もスコアが高い組($X_k,\Pi_{X_k}^l$)を見つける．このとき，最も高いスコアが改善していなければ終了．
$$k,l\leftarrow{\rm argmax}_ {i,j}{\rm Score}(X_i,\Pi_{X_i}^j)$$
$\Pi_{X_k}^l$を$\Pi_{X_k}$に置き換える，つまり，$\Pi_{X_k}\leftarrow\Pi_{X_k}^l$．その後，2に戻る．

BOAでは，評価指標Scoreとして，ベイズ情報量規準(Bayesian Information Criterion; BIC)やK2メトリックが用いられる(ことが多いと思う)．
以下には，近年の論文でよく用いられているBICの計算式を示す．

{\rm BIC}(X_i|\Pi_{X_i})=-H(X_i|\Pi_{X_i})\lambda - 2^{|\Pi_{X_i}|}\frac{\log_2\lambda}{2},\\\
H(X_i|\Pi_{X_i})=-\sum_{x_i,\pi_{x_i}}P(x_i,\pi_{x_i})\log_2P(x_i|\pi_{x_i}),

ここで，$\lambda$は集団サイズ，$H(X_i|\Pi_{X_i})$は$\Pi_{X_i}$のもとでの$X_i$の条件付きエントロピー，$x_i$は$X_i$の取りうる値(2値であれば0と1)，$\pi_{x_i}$は$\Pi_{X_i}$の取りうる値を意味する．

K2アルゴリズムによって各ノードの親ノードの集合を決定することができるので，実際にノード間にエッジを張ることでネットワーク構造が構築できる．

ノードの条件付き確率の推定

ベイジアンネットワークの各ノードの条件付き確率$P(x_i|\pi_{x_i})$は，最尤推定によって求める．
上記の3次元ビット列のベイジアンネットワークと下図(左)の集団が与えられたときの条件付き確率の推定結果は下図(右)のようになる．

まず，ルートノード$X_1$の確率を集団内の頻度情報をもとに決定する．
次に，$X_1$の子ノードである$X_2$について$X_1$の条件付き確率として推定する．
同様に葉ノード$X_3$も推定する．

ベイジアンネットワークからの個体の生成

個体は，構築したネットワーク構造$B$をルートノードから葉ノードまで順に辿り，推定した条件付き確率$P(X_i|\Pi_{X_i})$を用いて値を決定していくことで生成できる．

上記の例では，$X_1$の値を$P(X_1)$に基づいて決める．ここでは，$X_1=0$と仮定する．
次に，$X_2$の値を$P(X_2|X_1=0)$に基づいて決める．ここでは，$X_2=0$と仮定する．
最後に，$X_3$の値を$P(X_3|X_1=0,X_2=0)$に基づいて決める．ここでは，$X_3=1$と仮定する．
以上より，最終的に$X=(X_1,X_2,X_3)=(0,0,1)$という個体が生成される．

アルゴリズム

実装

プログラム全体はgithubに置いてある．

環境

Ubuntu18.04
Python 3.6.5
anaconda Command line client (version 1.7.2)
conda 4.8.3
numpy 1.18.1

githubの方にSingularityの環境ファイルを置いてあるので，それをビルドしても動く(動作確認済み 2019/12/19時点)．

BOAの実装を以下に示す．
基底クラスのEDABaseクラスに関しては分布推定アルゴリズムを学ぶ (2) ～PBIL～を参照

import numpy as np

from eda.optimizer.eda_base import EDABase
from eda.optimizer.metric import *
from eda.utils import estimate_cpt


class BOA(EDABase):
    def __init__(self, categories, selection, replacement,
                 lam=500, k=None, criterion="bic", theta_init=None):
        super(BOA, self).__init__(categories, lam=lam, theta_init=theta_init)
        assert k is None or 0 < k <= self.d
        assert criterion in ["aic", "bic", "k2"], \
            "Evaluation criterion takes one of three criterions, aic, bic, and k2."
        self.selection = selection
        self.replacement = replacement
        self.constraint_k = k
        if criterion == "aic":
            self.metric = AIC
        elif criterion == "bic":
            self.metric = BIC
        else:
            self.metric = K2

        self.population = None
        self.evaluation = None
        self.structure_mat = None
        self.sampling_order = None
        self.prob_mat = None

    def update(self, x, evals, range_restriction=False):
        x, evals = self._preprocess(x, evals)
        if self.population is None:
            self.population = x
            self.evaluation = evals
            self.lam = int(self.lam * self.replacement.replace_rate)
        else:
            self.population, self.evaluation = self.replacement(self.population,
                                                                self.evaluation,
                                                                x,
                                                                evals)
        x, evals = self.selection(self.population, self.evaluation)

        self.k2_algorithm(np.argmax(x, axis=2))


    def k2_algorithm(self, x):
        metric = self.metric(x, self.C)
        self.structure_mat = np.zeros((self.d, self.d))
        inheritance_mat = np.zeros((self.d, self.d))
        node_scores = np.zeros(self.d)
        child_counter = np.zeros(self.d)
        for node_idx in range(self.d):
            parents_idx = self.structure_mat[node_idx] != 0
            node_scores[node_idx] = metric.local_score(node_idx, parents_idx)

        min_gain = 0
        candidate_gains = np.zeros((self.d, self.d)) + min_gain
        for i in range(self.d):
            for j in range(i+1, self.d):
                parents_idx = self.structure_mat[i] != 0
                parents_idx[j] = True
                new_score = metric.local_score(i, parents_idx)
                candidate_gains[i, j] = new_score - node_scores[i]
                candidate_gains[j, i] = candidate_gains[i, j]
        # build graph structure with greedy algorithm
        while self.constraint_k is None or self.constraint_k > 0:
            child, parent = np.unravel_index(np.argmax(candidate_gains), candidate_gains.shape)
            if candidate_gains[child, parent] == min_gain:
                break
            self.structure_mat[child, parent] = 1
            node_scores[child] += candidate_gains[child, parent]

            candidate_gains[child, parent] = min_gain
            candidate_gains[parent, child] = min_gain

            inheritance_mat[child, parent] = 1

            new_ancestors = [parent]

            for i in range(self.d):
                if inheritance_mat[parent, i] == 1 and inheritance_mat[child, i] == 0:
                    new_ancestors.append(i)
                    inheritance_mat[child, i] = 1
                    candidate_gains[i, child] = min_gain

            descendents = np.where(inheritance_mat[:, child] == 1)[0].tolist()
            while len(descendents) > 0:
                node = descendents.pop(0)

                node_updated = False
                for ancestor in new_ancestors:
                    if inheritance_mat[node, ancestor] == 0:
                        inheritance_mat[node, ancestor] = 1
                        candidate_gains[ancestor, node] = min_gain
                        node_updated = True
                if node_updated:
                    descendents.extend(np.where(inheritance_mat[:, node] == 1)[0])
                    descendents = list(set(descendents))
            # update candidate gains for affected child node
            for i in range(self.d):
                if self.structure_mat[child, i] == 0 and inheritance_mat[i, child] == 0 and child != i:
                    parents_idx = self.structure_mat[child] != 0
                    parents_idx[i] = True
                    new_score = metric.local_score(child, parents_idx)
                    candidate_gains[child, i] = new_score - node_scores[child]

            child_counter[child] += 1
            if self.constraint_k is not None and child_counter[child] == self.constraint_k:
                candidate_gains[child] = min_gain

        self.prob_mat = []
        for i in range(self.d):
            parents_idx = self.structure_mat[i] == 1
            cpt = estimate_cpt(x[:, i], x[:, parents_idx], base=self.Cmax)
            self.prob_mat.append(cpt)
        # deceide sampling order of each node
        tmp_structure = self.structure_mat.copy()
        s = list(np.where(np.sum(tmp_structure, axis=1) == 0)[0])
        self.sampling_order = []
        while len(s) > 0:
            n = s.pop(0)
            self.sampling_order.append(n)

            candidates = np.where(tmp_structure[:, n] == 1)[0]
            for i in range(len(candidates)):
                m = candidates[i]
                tmp_structure[m, n] = 0
                if np.all(tmp_structure[m] == 0):
                    s.append(m)
        if np.any(tmp_structure != 0):
            print("error. graph has cycles")
            exit(-1)

        return metric.score(self.structure_mat)

    def sampling(self):
        # random sampling, only first generation
        if self.prob_mat is None:
            return super(BOA, self).sampling()
        # sample by using each probability of bayesian network
        else:
            c = np.zeros((self.d, self.Cmax), dtype=bool)
            for i in self.sampling_order:
                rand = np.random.rand()
                parents = self.structure_mat[i] == 1
                val = np.argmax(c[parents], axis=1)
                if len(val) == 0:
                    theta = self.prob_mat[i]
                else:
                    theta = self.prob_mat[i][tuple(val)]
                    if theta.shape[-1] == 0:
                        theta = theta.squeeze(-1)
                cum_theta = theta.cumsum()
                _c = (cum_theta - theta <= rand) & (rand < cum_theta)
                c[i, _c] = True
            return c

実験：BOAの性能評価

設定

目的関数はOne-Max問題，3-Deceptive問題を用いる．
One-Max問題は$D=100$，3-Deceptive問題は$D=30$とする．
また，3-Deceptive問題のユーザーパラメータは$d=0.1$とする．

その他のユーザーパラメータは以下の通り．

パラメータ名	値
集団サイズ$\lambda$	{100, 500, 1000, 1500}
選択手法	上位選択
選択割合	集団サイズの50%
置換方法	下位置換
置換割合	集団サイズの50%
BNの評価基準	BIC

以下のいずれかを満たした場合，最適化を終了する．

大域的最適解をサンプルできたとき(成功)
確率分布が収束したとき(失敗)
上限世代数($2\times 10^3$)に達したとき(失敗)

以上の設定で，乱数のシードを変更して10試行行う．

観察指標

以下の指標を記録する．

成功試行数(大域的最適解を見つけられた回数)
平均評価回数($\pm$標準偏差)
最適化終了時の最良評価値の平均($\pm$標準偏差)

ここで平均評価回数は，成功試行時のみのものとする．

結果

One-Max問題($D=100$)をBOAで最適化したときの結果を以下の表に示す．

集団サイズ	成功試行数	平均評価回数($\pm$標準偏差)	平均最良評価値($\pm$標準偏差)
100	0/10	-	94.60$\pm$2.06
500	10/10	5125.00$\pm$167.71	100.00$\pm$0.00
1000	10/10	9950.00$\pm$415.33	100.00$\pm$0.00
1500	10/10	15000.00$\pm$474.34	100.00$\pm$0.00

3-Deceptive問題($D=30$)をBOAで最適化したときの結果を以下の表に示す．

集団サイズ	成功試行数	平均評価回数($\pm$標準偏差)	平均最良評価値($\pm$標準偏差)
100	0/10	-	9.45$\pm$0.15
500	9/10	3444.44$\pm$421.27	9.98$\pm$0.06
1000	9/10	5722.22$\pm$415.74	9.99$\pm$0.03
1500	10/10	7875.00$\pm$768.52	10.00$\pm$0.00

いずれの問題もECGAの結果と似たような傾向になった．

One-Max問題では変数間の依存関係を考慮する必要がないので，どんなベイジアンネットワークでもある程度解けてしまい，集団サイズが500程度で最適化が可能という結果になった．

3-Deceptive問題では，安定して最適解を見つけるためにはECGA以上の集団サイズが必要となることがわかった(ECGAは500)．
ただ，同じ集団サイズで比較するとECGAよりも評価回数が少ないので，小さな集団サイズで安定して最適化できるようになれば，より優れる手法になるかもしれない．

考察

以下に，30次元の3-Deceptive問題を集団サイズ1500で最適化したときのベイジアンネットワークのネットワーク構造の一例を隣接行列によって示す．
縦軸と横軸の値が問題の次元に対応しており，色のついている部分が行変数と列変数の間に依存関係を検知している部分となる．
図より，対角線上に3変数ずつの塊で依存関係を検知しており，3-Deceptive問題の問題構造を把握できていることが確認できた．

まとめ

BOAはベイジアンネットワークによって問題構造を表現
騙し構造を持った問題に強く，変数間依存性を考慮しない手法では解けない問題を解くことができる
最適解を安定して発見するには，他手法と比べて大きな集団サイズを設定する必要がある
変数間に依存関係の存在しない問題での性能はPBILやCGAと比べて劣る