粒子群最適化と捕食者被食者最適化

# はじめに
この記事は[明治大学Advent Calendar 2018](https://qiita.com/advent-calendar/2018/meiji_univ)の7日目の記事です。7日目にして3度目の登場です。ただ、7日午前4時現在で18/25が埋まっていて、ユニーク投稿者数が9人となっているので、わくわくしてきました。

# 目標
**粒子群最適化と捕食者被食者の紹介と数値計算結果をもとに2手法の簡単な比較を行う。**

- **注意1** : 本記事にて使用したプログラムはProcessing(ほぼJava)で書かれています。
- **注意2** : 本記事で紹介する捕食者被食者最適化はあくまでも一つの例であり、ここでは一番シンプルなものを紹介している。他にも提案手法があるようなので、他が知りたければ[粒子群最適化や捕食者被食者最適化の亜種](https://qiita.com/ryunryunryun/private/e32837077e631c9e482c#%E7%B2%92%E5%AD%90%E7%BE%A4%E6%9C%80%E9%81%A9%E5%8C%96%E3%82%84%E6%8D%95%E9%A3%9F%E8%80%85%E8%A2%AB%E9%A3%9F%E8%80%85%E6%9C%80%E9%81%A9%E5%8C%96%E3%81%AE%E4%BA%9C%E7%A8%AE)を参照のこと。

# 粒子群最適化
## 概要
粒子群最適化(Particle Swarm Optimization; PSO)とはメタヒューリスティックな最適化手法の一つで、その名の通り粒子が群れを成して最適な点を探していく手法である。例えるなら動物の群れが、良さそうな餌場を見つけたリーダーに向かって行動をするような方法である。

## 粒子の更新則
まず、最適化した関数を$f: \mathbb{R^n} \to \mathbb{R}$とする。この$f$を最小にするような$x$を$\mathbb{R}^n$の中で見つけたい。そこで粒子群最適化では、最適な点を探索していくような粒子を$N$個用意する。そのために各粒子は情報として、その位置と速度を保持している。また、粒子の移動方法には以下の3つの特徴がある。

> **粒子の特徴**
1. 各粒子は全ての粒子の内、その時点で一番適した場所にいる粒子(リーダー)の居場所を知っていて、その粒子の方へと向かいたい。
2. 各粒子は記憶を持っていて、自分が通過してきた場所の中で一番適していた場所を覚えていて、その場所に行きたい。
3. 各粒子には慣性の法則が成り立っており、急に止まることはできない。

以上の特徴をもとに、各粒子の位置と速度の更新則を以下のようなものにする。

```math
\begin{align*}
v^{k + 1}_i & = c_1^k(x_i^k - x_{leader}^k) + c_2^k(x_i^k -\bar{x}_{i}^k) + w_i^kv_i^k \\
x^{k + 1}_i & = x_i^k + v_i^{k + 1}
\end{align*}
```

ここで、$x_i^k$と$v_i^k$はそれぞれ$i$番目の粒子の$k$ステップ目の位置と速度を表している。速度の更新則に出てくる各パラメータは、1項目は$c_1^k$があるパラメータ$c_1$に対して$0 \sim c_1$の間の乱数、$x_{leader}^k$は$k$ステップ目において一番適した場所にいる粒子の場所、すなわち

```math
x_{leader}^k \in \underset{x_i^k}{\operatorname{argmin}} \{f(x_i^k) \mid i = 0, \ldots, N - 1\}
```

である。2項目は、$c_2^k$があるパラメータ$c_2$に対して$0 \sim c_2$の間の乱数、$\bar{x}_i^k$が$i$番目の粒子の$k$ステップまでで一番適していた場所、すなわち

```math
\bar{x}_i^k \in \underset{x_i^j}{\operatorname{argmin}} \{f(x_i^j) \mid j = 0, \ldots, k\}
```
である。3項目の$w$はあるパラメータ$w$に対して$0 \sim w$の間の乱数で慣性の度合いを表す。

**それぞれの項がそれぞれの移動方法の特徴を表している。**

## アルゴリズム
粒子群最適化のアルゴリズムを記す。

1. $N$個の粒子を、初期位置をランダムに$x_0^0, \ldots, x_{N - 1}^0$、初速度もランダムに$v_0^0, \ldots, v_{N - 1}^0$として用意する。$(k = 0)$
2. 各粒子の適応度$f(x_i^k), \ (i = 0, \ldots, N - 1)$を求めて$x_{leader}^k$と$\bar{x}_i^k$を求める。
3. 上記の更新則をもとに$v_0^{k + 1}, \ldots, v_{N - 1}^{k + 1}$と$x_0^{k + 1}, \ldots, x_{N - 1}^{k + 1}$を求める。
4. $\left|x_{leader}^k - x_{leader}^{k + 1}\right|$が十分小さくなるまで、2と3を繰り返す。

# 捕食者被食者最適化
## 概要
捕食者被食者最適化とは粒子群最適化に対して粒子(被食者)を追いかけていくような捕食者粒子を付け足した最適化手法である。あとで数値計算例を見るが、粒子群最適化では局所最適解にすら到達しない場合があるのに対して、捕食者被食者最適化では少なくとも局所最適解には到達しやすいようになっている。

## 粒子の更新則
被食者粒子の速度の更新則は粒子群最適化の3つの特徴に加えて、捕食者粒子が近づいてきたらパニックに陥るという特徴がある。また、捕食者粒子はリーダー被食者粒子へまっしぐらに移動していく。

$x_i^k$と$v_j^k$を$k$ステップ目での$i$番目の被食者粒子の位置と速度として、$x_{pred}^k$と$v_{pred}^k$を捕食者粒子の位置と速度とする。このとき、更新則は以下のように表される。

```math
\begin{align*}
v^{k + 1}_i & = c_1^k(x_i^k - x_{leader}^k) + c_2^k(x_i^k -\bar{x}_{i}^k) + w_i^kv_i^k + D(|x_i^k - x_{pred}^k|)K_i^k\\
x^{k + 1}_i & = x_i^k + v_i^{k + 1} \\
v_{pred}^{k + 1} & = c_{pred}^k (x_{pred}^k - x_{leader}^k) \\
x^{k + 1}_{pred} & = x_{pred}^{k} + v_{pred}^k
\end{align*}
```

各種パラメータの内、被食者粒子の速度の3項目までは粒子群最適化のものと同じである。4項目に関しては$K_i^k$は$K$をパラメータとして半径$K$の球面状のランダムな点で、$D$は$b$を正のパラメータとして以下のように定義される関数である。

```math
D(|x_i^k - x_{pred}^k|) = e^{-b |x_i^k - x_{pred}^k|}
``` 
であり、捕食者が近づいてきたらパニックに陥いる様子を表している。また、捕食者の速度に関しては、1項目の$c_{pred}^k$はパラメータ$c_{pred}$に対して$0 \sim c_{pred}$の間の乱数である。

## アルゴリズム

捕食者被食者最適化のアルゴリズムを記す。

1. $N$個の被食者粒子と1つの捕食者粒子を、初期位置をランダムに$x_0^0, \ldots, x_{N - 1}^0, x_{pred}^0$、初速度もランダムに$v_0^0, \ldots, v_{N - 1}^0, v_{pred}^0$として用意する。$(k = 0)$
2. 各被食者粒子の適応度$f(x_i^k), \ (i = 0, \ldots, N - 1)$を求めて$x_{leader}^k$と$\bar{x}^k_i$を求める。
3. 上記の更新則をもとに$v_0^{k + 1}, \ldots, v_{N - 1}^{k + 1}, v_{pred}^k$と$x_0^{k + 1}, \ldots, x_{N - 1}^{k + 1}, x_{pred}^k$を求める。
4. 所望の結果が得られるまで、2と3を繰り返す。

# 数値計算例
プログラムは現在ブラッシュアップ中です。その内Githubにでも上げて、リンクをこちらに貼ります。 2手法に対してパラメータは以下のように設定した。

```math
c_1 = 0.1, \ c_2 = 0.3, \ w = 0.1 \\
\ K = 0.5, \ b = 10, \ c_{pred} = 0.1
```

また、目的関数は次のような$\mathbb{R}^2$上のRosembrock関数というテストにしばしば用いられる関数を用いた。なお、最適解は$x_{opt} = (1, 1)$である。$$f(x, y) = 100 (y - x)^2 + (1 - x)^2, \quad (-2.048 \le x, y \le 2.048).$$


## 粒子群最適化
**ピンクの三角が最適解(1, 1)であり、赤い丸が各粒子である。また、色が濃いところほど$f$の値が低い。**

- 粒子数10の場合: 最適化ではない場所で静止してしまった。
![out1.gif](https://qiita-image-store.s3.amazonaws.com/0/217692/017e6c1f-9aca-e14e-22e5-e2899901e8fb.gif)

- 粒子数1000の場合: 数の暴力で最適化に行った。
![out2.gif](https://qiita-image-store.s3.amazonaws.com/0/217692/1fd0703c-9c71-e5c4-f105-49b03ee01f0e.gif)


## 捕食者被食者最適化

- **青が捕食者粒子。**　被食者粒子数は10: 最初は最適解とは程遠い場所に集まったが、捕食者に追いやられて最適解に近づいた。しかし、パニックは収まらず、最適解の近くをわらわらしている。
![out3.gif](https://qiita-image-store.s3.amazonaws.com/0/217692/23f0fb33-d257-7751-4025-777f310fdbc1.gif)


# 2手法の比較
粒子群最適化では粒子数がものを言う。この手法では、粒子がリーダー粒子に近づく中で最適な点を見つけられるかどうかが、大域的な最適解の探索や局所最適解から抜け出す可能性に繋がる。そのため、粒子数が少ないと、それだけ探索範囲が狭まり、局所的な最適解にすら行き着かない場合が出てくる。しかし、粒子が増えればそれだけ良いリーダーが生まれやすく、大域的な最適解に行き着く確率が上がる。

捕食者被食者最適化では粒子数が最適解の探索にそこまで影響しない。捕食者粒子が近くにいない場合、被食者粒子達はほぼ粒子群最適化と同様の振る舞いを示し、リーダーのもとへと集まろうとする。このとき、リーダーは最適解の近くにいるとは限らず、全く見当違いの場所にいるかもしれない。しかし、捕食者粒子が近づいてくると、皆パニックを起こし、ほぼランダムサーチのような状態になる。こうすることで、見当違いの場所で落ち着いていた可能性のある群れが、より(局所)最適解の方向へ追いやられていくのである。

# 粒子群最適化や捕食者被食者最適化の亜種
- 白髪 一将, 原 章, 高濱 徹行, 「捕食者被食者の関係を導入したHeterogeneous Particle Swarm Optimization」, 日本知能情報ファジィ学会 ファジィ システム シンポジウム 講演論文集, 2012, 28 巻, 第28回ファジィシステムシンポジウム, p.727-732, https://www.jstage.jst.go.jp/article/fss/28/0/28_727/_article/-char/ja/
    - 今回紹介した捕食者被食者最適化が載っている。他にもHeterogeneous Particle Swarm Optimizationという、粒子毎の更新則に個性をつけた粒子群最適化の亜種も紹介されており、これを元にした捕食者被食者最適化を提案している。

- 犬飼規雄, 井上拓也, 上手洋子, 西尾芳文, 「捕食者を加えた粒子群最適化での捕食範囲の調査」, 電子情報通信学会 信学技報, vol. 113, no. 271, NLP2013-90, pp. 109-112, 2013年10月. https://www.ieice.org/ken/paper/20131028yB6k/
    - 本記事の捕食者被食者最適化は「捕食者被食者」と言いながらも追いかけ追いかけられの関係のみで捕食被食の様子はなかった。しかし、この論文では捕食者粒子が実際に被食者粒子を食べる。食べられた粒子は消え、その分またランダムな位置に復活する。

# 結論

粒子群最適化と捕食者被食者最適化を紹介し、その性能について簡単な比較を行った。