分布推定アルゴリズムを学ぶ (3) ～UMDA～

※この記事は以下の記事を読んでいることを想定しています．
[分布推定アルゴリズムを学ぶ (1) ～概要～](https://qiita.com/_tetsuya_/items/0cc2123f87007e91bd1c)
[分布推定アルゴリズムを学ぶ (2) ～PBIL～](https://qiita.com/_tetsuya_/items/1612e616964091d46787)


## Univariate Marginal Distribution Algorithm(UMDA)とは
変数間の依存性を考慮しない離散最適化手法の一つ．
主にビット列を対象としている．

分布推定アルゴリズム(Estimation of Distribution Algorithms; EDAs)に属する．


## 概要
UMDAではEDAsの基本的な流れに沿って最適化を進め，確率モデルM<sup>t</sup>部分に単変量周辺分布と呼ばれる分布を用いてモデルの推定を行う．

個体$\boldsymbol{x}$を$\boldsymbol{x}=(x_1,x_2,...,x_D)\in \\{0,1\\}^D$のD次元ビット列と仮定する．
有望な解集団S<sup>t</sup>内の各次元での0と1の出現頻度を新たな個体を生成するための確率分布M<sup>t</sup>として用いる．
そのため，$p(x_i)$を$x_i$が1になる確率とすると，新たな個体$\boldsymbol{x}$を生成する確率分布は以下で表現される．

```math
p(\boldsymbol{x})=\prod_{i=1}^{D}p(x_i)
```

UMDAは線形問題に対して完全であることが示されている(らしい)．

## 実装
プログラム全体は[github](https://github.com/e5120/EDAs)に置いてある．

### 環境
- Ubuntu 18.04.4 LTS
- Python 3.7.4
- anaconda Command line client (version 1.7.2)
- conda 4.8.1
- numpy 1.17.2

[github](https://github.com/e5120/EDAs)の方に[Singularity](https://sylabs.io/guides/3.2/user-guide/index.html)の環境ファイルを置いてあるので，それをビルドしても動く(動作確認済み 2019/12/19時点)．

UMDAの実装を以下に示す．
基底クラスの`EDABase`クラスに関しては[分布推定アルゴリズムを学ぶ (2) ～PBIL～](https://qiita.com/_tetsuya_/items/1612e616964091d46787)を参照．

`selection`変数は任意の選択方式(トーナメント選択やルーレット選択等)のオブジェクトを指す．
`c_one`変数はone-hot形式の個体の集団を想定しており，`(集団サイズ,次元数,one-hot)`の形を取る．
`fxc`変数は`c_one`変数の各個体に対応する評価値のリスト．

```python
import numpy as np

from optimizer.eda_base import EDABase


class UMDA(EDABase):
    def __init__(self, categories, lr, selection, lam=64, theta_init=None):
        super(UMDA, self).__init__(categories, lam=lam, theta_init=theta_init)
        # only binary strings
        assert self.Cmax == 2
        assert 0 < lr < 1
        self.selection = selection
        self.lr = lr

    def update(self, c_one, fxc, range_restriction=False):
        self.eval_count += c_one.shape[0]
        # sort by fitness
        idx = np.argsort(fxc)
        # store best individual and evaluation value
        if self.best_eval > fxc[idx[0]]:
            self.best_eval = fxc[idx[0]]
            self.best_indiv = c_one[idx[0]]
        # apply selection
        if self.selection is not None:
            c_one, fxc = self.selection(c_one, fxc)
        # update probability vector
        self.theta[:, -1] += self.lr * (np.mean(c_one[:, :, -1], axis=0) - self.theta[:, -1])
        self.theta[:, 0] = 1 - self.theta[:, -1]
        # if range_restriction is True, clipping
        for i in range(self.d):
            ci = self.C[i]
            theta_min = 1.0 / (self.valid_d * (ci - 1)) if range_restriction and ci > 1 else 0.0
            self.theta[i, :ci] = np.maximum(self.theta[i, :ci], theta_min)
            theta_sum = self.theta[i, :ci].sum()
            tmp = theta_sum - theta_min * ci
            self.theta[i, :ci] -= (theta_sum - 1.0) * (self.theta[i, :ci] - theta_min) / tmp
            self.theta[i, :ci] /= self.theta[i, :ci].sum()
```

## 実験：UMDAの性能評価
### 設定
目的関数はOne-Max問題を用い，$D=100$とする．

その他のユーザーパラメータは以下の組み合わせ全てを試す．

|パラメータ名 |値 |
|---|---|
|集団サイズ：$N$ |{8, 32, 64} |
|学習率：$lr$ |{0.1, 0.3, 0.5} |
|クリッピング処理 | {False, True} |
|選択方式|トーナメント選択(トーナメントサイズ=2) |

以下のいずれかを満たした場合，最適化を終了する．
1. 大域的最適解をサンプルできたとき(成功)
2. 確率分布が収束したとき(失敗)
3. 上限世代数($2\times 10^3$)に達したとき(失敗)

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

### 観察指標
以下の指標を記録する．

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

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

```math
評価回数 = 集団サイズ \times 世代数
```

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

|集団サイズ |学習率$lr$| クリッピング処理| 成功試行数 | 平均評価回数($\pm$標準偏差)|平均最良評価値($\pm$標準偏差)|
|---|---|---|---|---|---|
|8 |0.1|False| 9/10 |1855.11$\pm$189.44|99.9$\pm$0.30|
|32 |0.1|False|10/10 | 6316.80$\pm$283.41|100.0$\pm$0.0|
|64 |0.1|False|10/10 |12467.2$\pm$441.37|100.0$\pm$0.0|
|8 |0.3|False| 0/10 |-|89.50$\pm$1.63|
|32 |0.3|False|9/10 |2385.78$\pm$113.0|99.9$\pm$0.30|
|64 |0.3|False|10/10 |4371.20$\pm$198.40|100.0$\pm$0.0|
|8 |0.5|False| 0/10 |-|81.00$\pm$3.10|
|32 |0.5|False|5/10 |1529.60$\pm$84.42|99.1$\pm$1.04|
|64 |0.5|False|10/10 |2848.00$\pm$128.80|100.0$\pm$0.0|
|8 |0.1|True| 10/10 |1983.20$\pm$152.77|100.0$\pm$0.0|
|32 |0.1|True|10/10 |6396.80$\pm$277.29|100.0$\pm$0.0|
|64 |0.1|True|10/10 |12531.20$\pm$429.13|100.0$\pm$0.0|
|8 |0.3|True| 10/10 |2205.60$\pm$490.89|100.0$\pm$0.0|
|32 |0.3|True|10/10 |2409.60$\pm$165.69|100.0$\pm$0.0|
|64 |0.3|True|10/10 |4467.20$\pm$180.57|100.0$\pm$0.0|
|8 |0.5|True| 10/10 |2237.60$\pm$590.45|100.0$\pm$0.0|
|32 |0.5|True|10/10 |1987.20$\pm$234.91|100.0$\pm$0.0|
|64 |0.5|True|10/10 |2963.20$\pm$218.07|100.0$\pm$0.0|

### 考察
集団サイズが小さくて学習率が大きかったり，集団サイズが大きくて学習率が小さいと探索がうまく行かなかったり非効率になることが確認できた．
パラメータの設定に評価回数が大きく影響されるため，パラメータの探索がある程度必要な手法であると考えられる．

集団サイズが8で学習率が0.3や0.5だと確率分布が誤った方向に収束してしまうため，成功試行数が0になっている．
しかし，早期の収束を回避するためにクリッピング処理を行うと全試行成功するようになり，クリッピング処理の有用性が示された．

## まとめ
1. 変数間の依存性を持たない問題は解ける
2. ユーザーパラメータの設定に対して性能が敏感に変化
3. PBILと比較すると同程度または少し劣るので，PBILで十分という感想

## 参考文献
1. [Marginal Distributions in Evolutionary Algorithms](https://www.researchgate.net/publication/2421623_Marginal_Distributions_in_Evolutionary_Algorithms)