粒子群最適化（Particle Swarm Optimization: PSO）をGoで実装してみた

粒子群最適化とは、群知能による最適化手法の一種です。この手のバイオミメティクス（生物模倣）によるアプローチは、私の学生時代の専門である材料工学でも非常に盛んに研究されていましたが、データサイエンスでも応用されているのを知って興味を持ったので、実装しながら勉強しました。

粒子群最適化（Particle Swarm Optimazation: PSO）

粒子群最適化とは、鳥や魚などの群れに見られる社会的行動のシミュレーションを基にモデル化されたヒューリスティックな最適解探索アルゴリズムです。

全体的な最適化の流れを説明します。まず最初に探索空間内で、決められた個数の粒子を全て初期化します。初期化とは空間内の位置をランダムに決めることを指します。

それぞれの粒子は、目的関数により適合度（どれだけその解が優れているか）を算出できます。全ての粒子をある移動ルールを元に動かしながら、適合度が高い位置を探していくのがPSOの探索方法です。

ここでPSOにおいて重要な移動ルールについて説明します。

粒子 $P_{i}$ の時刻 $t$ から $t+1$ の移動速度 $\vec{v}_{i}(t+1)$ は以下の式で表されます。

f:id:hirotsuru314:20190409205526p:plain:w500

ここで重要なのが、グローバルベスト $\vec{g}(t)$ とパーソナルベスト $\vec{p}_{i}(t)$ です。
グローバルベストとは全ての粒子群で過去の探索も含めて最も適合度が高い位置、パーソナルベストとは特定の1粒子におけるこれまでの探索の中で最も適合度が高い位置を示します。これは繰り返しの探索の中で随時更新されていく変数です。
$I$ は慣性係数と呼ばれる定数で、前の移動速度がどれだけ維持されるかを表します。 $A_{g}$ および $A_{p}$ は加速係数と呼ばれる定数です。

この数式からわかるように粒子群の中ではグローバルベストの情報が共有されていて、その方向に向かうように移動速度を調整しながら、より良い解の探索を繰り返します。これはあたかも餌を見つけた鳥が仲間にその情報を伝えて、仲間が餌のある位置に集まってくるようです。

作ったもの

github.com

『Evolutionary Algorithm implemented in Go』の略で名前を付けました。進化計算全般に興味があるので、ここに実装しながら勉強したいと思っています。

最適化問題を解いてみる

今回は下のシンプルな目的関数の最小値を最適解として、PSOでパラメータを探索してみました。

f:id:hirotsuru314:20190409212421p:plain:w120

$N=2$ としたときの関数のグラフを描くと以下のようになります。

f:id:hirotsuru314:20190409215708p:plain

左側が3Dに描いたもので、右側が2次元に写像したもの（左のグラフを上から見たもの）です。
グラフや数式から直感的にわかるように、この関数は $(x_1,x_2)=(0,0)$ で最小値 $0$ を取るため、今回の探索の目的は、 $(x_1,x_2)=(0,0)$ を見つけることになります。

サンプルコードは以下の通りです。

package main

import (
    "github.com/tsurubee/eago"
    "log"
)

func objectiveFunc(x []float64) float64 {
    return x[0] * x[0] + x[1] * x[1]
}

func main() {
    pso := eago.NewDefaultPSO()
    pso.NParticle =  5
    pso.NStep = 20
    pso.Min = -20
    pso.Max = 10

    if err := pso.Minimize(objectiveFunc, 2); err != nil {
        log.Fatal(err)
    }
}

上のコードをmain.goとして以下のように実行します。

$ go run main.go
Step   0: Fitness=8.124 Position=[-2.168 1.851]
Step   1: Fitness=8.124 Position=[-2.168 1.851]
Step   2: Fitness=7.394 Position=[-1.618 2.186]
Step   3: Fitness=6.958 Position=[1.068 2.412]
Step   4: Fitness=6.958 Position=[1.068 2.412]
Step   5: Fitness=4.809 Position=[-2.106 0.611]
Step   6: Fitness=2.215 Position=[-1.052 1.053]
Step   7: Fitness=1.535 Position=[0.580 1.095]
Step   8: Fitness=0.163 Position=[0.023 0.403]
Step   9: Fitness=0.163 Position=[0.023 0.403]
Step  10: Fitness=0.055 Position=[-0.175 0.156]
Step  11: Fitness=0.055 Position=[-0.175 0.156]
Step  12: Fitness=0.055 Position=[-0.175 0.156]
Step  13: Fitness=0.055 Position=[-0.175 0.156]
Step  14: Fitness=0.055 Position=[-0.175 0.156]
Step  15: Fitness=0.055 Position=[-0.175 0.156]
Step  16: Fitness=0.055 Position=[-0.180 0.151]
Step  17: Fitness=0.044 Position=[0.019 0.210]
Step  18: Fitness=0.013 Position=[-0.111 0.023]
Step  19: Fitness=0.008 Position=[-0.083 -0.031]

出力結果を見てみると、Stepが進むにつれて、Fitnessが $0$ 、パラメータの解が $(x_1,x_2)=(0,0)$ に収束しようとしていることがわかります。
粒子が速度を調整しながら最適解に収束して行く様子は下のサイトのアニメーションが非常にわかりやすかったです。

Particle swarm optimization - Wikiwand