PySwarms による粒子群最適化 - ご注文はリード化合物ですか？〜医薬化学録にわ〜

世の中はワクチンや~~東京オリンピックや~~ AlphaFold2 で盛り上がっていますが、久々に書く記事は流行とは少し離れた内容です。

最適化手法というと、基本的な最急降下法や、それこそ流行りのベイズ最適化などがありますが、その一つに、粒子群最適化（partial swarm optimization: PSO）というものがあります。PSO は群知能（例えば渡り鳥一匹一匹はアホでも、集団で動くことで最適な行動を選択できることに相当する。~~コロナで最適な行動を取れない人類とは違う~~）の一種であり、1995 年に Kennedy と Eberhart により報告されています。
PSO では、グローバルな最適解を見出すのに、たくさんの点を用います（粒子群と呼ばれる所以です）。
各々の点は、普通の最急降下法のように、局所最適な方向に動こうとします。しかし、座標の更新にあたっては、局所的に最適な方向だけでなく、全ての点の中で最も良い座標も考慮します。そのため、PSO の最適化関数は、慣性項、ローカル項、グローバル項の 3 つの項から成り立っています。慣性項は、座標を更新（移動）する基本的なスピードを、ローカル項は、局所最適解を、グローバル項は、粒子群全体における最適値を考慮する役割をそれぞれ持っています（ざっくりとした理解なので、厳密な理解は原著論文をあたることをお勧めします）。

PSO ですが、Python で PySwarms というライブラリがあり、pip で簡単に入れることができます。

pyswarms.readthedocs.io

pip install pyswarms

PySwarm による PSO は動かすだけならかなり簡単に使えます。

import numpy as np
import pandas as pd
import pyswarms as ps

# 最適化したい関数を定義（10, -3, 3 で最小値 0 になる）
def object_function(x):
    y = (x[:, 0] - 10) ** 2 + (x[:, 1] + 3) ** 2 + (x[:, 1] + x[:, 2]) ** 2
    return y

# 最適化関数の重み(hyperparameter)を決める
options = {"c1": 0.5, "c2": 0.3, "w":0.9}

# 乱数固定
np.random.seed(0)

# PSO を実行するインスタンスを設定する
optimizer = ps.single.GlobalBestPSO(n_particles = 30, dimensions = 3, options = options)

# Perform optimization
cost, pos = optimizer.optimize(object_function, iters = 100)

# 最適化された値を表示
print(cost)

# 最適化された座標を表示
print(pos)

今回は 3 次元座標を入力すると数字が返ってくる関数を定義し、その最小値を探索することを考えます。
options で 3 つのパラメータを決めていますが、c1 がローカル項の重み、c2 がグローバル項の重み、w が慣性項の重みをそれぞれ示しています。w は 1 より大きくすると挙動がおかしくなるので、0~1 の範囲にしましょう。細かい範囲を探索するときは、値をオーダーレベルで小さくした方が良い気がします。c1 と c2 も重要なはずなのですが、筆者が変更した範囲では、あまり結果は変わりませんでした。

PySwarms は numpy に対応しているため、np.random.seed で乱数を固定し、結果を再現することができます（公式ドキュメントには書いていなかったが、明記しておいて欲しい）。

PySwarms では GlobalBestPSO と LocalBestPSO の二つのメソッドが使えますが、筆者が使ってみた限りでは差がありませんでした。
インスタンスを定義する際に、PSO の条件を色々設定できます。確実に使うのは、粒子数（n_particles）と dimension, option です。粒子数は、計算資源が許す限り、大きくしておいて損はないと思います。計算時間も、線形増加はしないです。dimension は、最適化したい目的関数のインプットの次元数です。今回は三次元データをインプットしているので、3 にしています。option では、PSO の最適化関数に関する設定を辞書として入力します。

optimizer.optimize で PSO を実行します。第一引数に最適化したい関数を入力します。関数は自分で定義できるので、機械学習モデルなどを用いることも可能です。また、値は小さくなるように更新されていくので、最大値を得たい場合は、関数にマイナスをつける必要があります。iters は繰り返しの回数で、これも単純に大きい方が良いですが、時間は線形増加します。

戻り値として、最適値（cost）と最適座標（pos）が帰ってきます。いずれも、最新の計算ではなく、今までの全ての計算の中で最適なものを算出しているので、計算のし過ぎによる悪影響はありません。

optimizer.cost_history で、最適値の履歴を時系列順のリストとして確認することができます。最適化が進んでいるか、プラトーに達しているかを判断できます。
optimizer.pos_history で、これまで計算した全ての座標を確認できます。時系列順のリストの中に、座標の numpy データが入っています。

状況によっては、初期座標を指定して計算したい場合もあると思います。座標を numpy 形式で用意し、optimizer を定義するときに init_pos で指定することができます。

# 初期座標を用意
initial_position = np.tile(np.arange(-15, 15) * 0.02 + 7, (3, 1)).T

# init_pos で初期座標を設定
optimizer = ps.single.GlobalBestPSO(n_particles = 30, dimensions = 3, options = options, init_pos = initial_position)

探索範囲を指定したい場合は、optimizer の bounds で範囲を指定できます。最小値と最大値のタプルを用意します。タプルの次元数は、探索次元に合わせる必要があります。

# 最適化関数の重み(hyperparameter)を決める
options = {"c1": 0.5, "c2": 0.3, "w":0.9}

# 乱数固定
np.random.seed(0)

# 探索範囲を指定。一つ目のタプルは最小値を、二つ目のタプルは最大値を示す。
bounds = ((7, -5, 4), (9, 0, 6))

# init_pos で初期座標を設定
optimizer = ps.single.GlobalBestPSO(n_particles = 30, dimensions = 3, options = options, bounds = bounds)

# Perform optimization
cost, pos = optimizer.optimize(object_function, iters = 100)

# 最適化された値を表示
print(pos)

# 最適化された座標を表示
print(cost)

最適化関数の重みである option は hyperparemeter ですが、gridsearch による最適化を PySwarms で行うことができます。

from pyswarms.utils.search import GridSearch

# Grid search による最適化関数各項の重みの最適化
options = {"c1": [0.3, 0.5, 0.8], "c2": [0.3, 0.5, 0.8], "w": [0.2, 0.3, 0.5]}
g_search = GridSearch(ps.single.GlobalBestPSO, n_particles = 10, dimensions = 3,
                   options = options, objective_func = object_function, iters = 20)

best_score, best_options = g_search.search()

PySwarms の基本的な使い方を今回はまとめてみました。PSO は現在でも時々使われている最適化手法なので、選択肢の一つに加えておくのも良いのではないかと思います。