Emileの備忘録

色々な事をだらだらと

ブラックボックス(ベイズ)最適化とハイパラの自動調整

この記事は学ロボアドベントカレンダー2023最終日の記事です。
学ロボアドカレ2023、多数のご参加ありがとうございました!!



さて、2023年も年の瀬が迫って参りました。
皆様いかがお過ごしでしょうか?
私はこの年末が学部生活の冬という現実に戦々恐々としています。

adventar.org




この記事では、ハイパーパラメータの自動調整に用いられるblackbox最適化を概説します。
blackbox最適化の中でも特にベイズ最適化を中心に解説し、optunaを用いた実装についても解説します。


環境としては、python 3.11.6を想定しています。




==============目次==============

===============================




1.ブラックボックス最適化とは

式で陽に表現できない関数をblackbox関数と言い、
blackbox関数に対する最適化問題blackbox最適化と言います。

以降、blackbox関数を f: \mathcal{X} \rightarrow \mathbb{R}と表す事にします。


特にblackbox最適化では以下の様な状況を考える事が多いです。

 \nabla f(x), \nabla^2 f(x)が観測できない
 f(x)の評価コストが高い
 f(x)の観測にノイズが乗る ( y = f(x)+\epsilonが観測される)


勾配やヘッセ行列などが得られないので、一般的な最適化手法を適用する事は出来ません。
blackbox最適化では、この様な関数 fに対して f(x)最小にする x \in \mathcal{X}を見つける事を考えます。



(例) blackbox最適化の例としては、何らかの問題のハイパーパラメータの自動調整が挙げられます。
この記事の後半では、ハイパラの自動調整に広く用いられるoptunaで実際にblackbox最適化を実装します。





1-1. Blackbox最適化の手法

主なBlackbox最適化の手法として、以下が挙げられます。

・ランダムサーチ/グリッドサーチ
文字通り、予め決めたルール/事前情報で総当りします。
これで問題無いなら一番手っ取り早いです。


ベイズ最適化
(この記事で主に解説する内容です。)


・進化計算
入力 x\in \mathcal{X}を個体とした進化計算によってblackbox最適化を行います。
GAやCMA-ESといった手法がよく使われます。

CMA-ESについてはこちらの記事の解説が分かり良いです。

horomary.hatenablog.com






2. ベイズ最適化

この記事では、特にblackbox最適化の中でも広く用いられているベイズ最適化を扱います。
ベイズ最適化では、blackbox関数 fに統計モデルを仮定し、 fを分布からのサンプルとして捉えます。


具体的には、以下の手続きで fを最小化する点 xを求めます。

1. 観測データ \mathcal{D}_n := \{(x_i, y_i)\}_{i=1}^nから統計モデルを更新
2. 統計モデルから獲得関数を求め、獲得関数を最大化する点 x'を次の入力として選択
3. 選択した点 x'を入力して、評価値 y'を観測し、 \mathcal{D}_{n+1} \leftarrow \mathcal{D}_n \cup \{(x', y')\}とする。


統計モデルは入力に対する fの分布の近似に使用し、
獲得関数は探索と活用のトレードオフを考慮した入力を定めるのに使用します。


強化学習などでもおなじみ探索と活用のトレードオフはここでも登場します。
勿論UCB系のアプローチも使われます。(今回はminimizeなのでLCB)






3. ベイズ最適化の手法

ベイズ最適化には、統計モデルと獲得関数の組み合わせで様々な種類の手法が存在します。
ここでは、ベイズ最適化でよく使われる手法を統計モデル/獲得関数のそれぞれで列挙します。


3-1. 統計モデル

まず統計モデルから概説します。
統計モデルを用いてblackbox関数 fをモデル化します。

ガウス過程回帰(GPR)

ガウス過程(GP)というのは、入力空間上のランダムな関数を定める確率過程です。

出力の平均値を定める関数  m_n: \mathcal{X} \rightarrow \mathbb{R}カーネル関数  k_n: \mathcal{X}\times \mathcal{X} \rightarrow \mathbb{R}から構成され、
入力系列  x_1,...,x_n \in \mathcal{X}に対する、関数 g: \mathcal{X} \rightarrow \mathbb{R}の出力系列が
 \mathbf{g}:= (g(x_1),.., g(x_n)) \sim \mathcal{N}(\mathbf{m}, \mathbf{K})となる時、関数 gガウス過程に従うと言います。
ただし、 \mathbf{m}_i = m_n(x_i) \mathbf{K}_{ij} = k_n(x_i, x_j)とします。


以降では、 \mathcal{D}_nが得られている時にGPの平均を表す関数を m_n(x)、分散を表す関数を k_n(x, x)として、
関数 gガウス過程に従う事を g \sim \mathcal{GP}(m_n, k_n)と表す事にします。

更に、観測 yに乗るノイズの分散を \sigma^2として、
 \sigma_n^2(x) := k_n(x, x) + \sigma^2を定義します。( y \sim \mathcal{N}(m_n(x), \sigma^2_n(x)))


ガウス過程回帰(GPR)は、ガウス過程と正規分布の再生性・同時分布の性質を用いて、
ガウス過程の m_n, k_nをオンラインで更新する事で関数回帰を行う手法です。

GPRでは、 fの事前分布を \mathcal{GP}(m_0, k_0)とし、観測ノイズ \epsilon \sim \mathcal{N}(0, \sigma^2)とすると、
 \mathcal{D}_n = \{(x_i, y_i)\}_{i=1}^nが観測された場合の事後分布は以下の様に表されます。


 f \sim \mathcal{GP}(m_n, k_n)
 m_n(x) = m_0(x) + \mathbf{k}_n(x)^\top (\mathbf{K}_n + \sigma^2 I_n)^{-1} ( \mathbf{y}_n - m_0(x))
 k_n(x, x') = k_0(x, x') - \mathbf{k}_n(x)^\top (\mathbf{K}_n + \sigma^2 I)^{-1} \mathbf{k}_n(x)

ただし、 \mathbf{k}_n (x) = (k_n(x_1, x),..., k_n(x_n, x))^\topとし (\mathbf{K}_{n})_{ij} = k_n(x_i, x_j)としています。


上記の様にガウス過程としてモデル化する事で、回帰する対象の分散(不確実さ)も表現出来るのがメリットです。
ただし、観測した点数 nについて更新にかかる計算量が O(n^3)と計算量が大きいです。




TPE(Tree-structured Parzen Estimator)

TPEでは、 P(x|y, \mathcal{D}_n)のみをモデル化して P(y|x, \mathcal{D}_n)を求めます。
獲得関数にはEI(詳細は後述)を用います。


具体的には、
 P(x|y, \mathcal{D}_n) := \left\{ \begin{array}{lll} l(x | \mathcal{D}_n) & (y < y^*) \\ g(x|\mathcal{D}_n) & (y \geq y^*)  \end{array} \right.としてモデル化し、
関数 l(x), g(x) \mathcal{D}_nを用いたParzen推定(カーネル密度推定法の一種)で求めます。
 P(y)は陽にモデル化しない事に注意して下さい。


ただし、 y^* P(y < y^*) = \gammaとなる様に与えた閾値とします。( \mathcal{D}_nの上位 \gammaとなる yとか)
optunaでは、 \gamma = 0.25がデフォルト値として採用されています。

 y y^*よりも小さいか否かを考えたいので、この様に分割しています。


このモデルを用いたEIを計算すると以下の様になります。
(TPEは獲得関数までセットなので、ここで解説します)

 \alpha_{EI}(x) = \int_{-\infty}^{y^*} (y^* -y)P(y | x,\mathcal{D}_n) dy
     = \int_{-\infty}^{y^*} (y^* - y)\frac{l(x, \mathcal{D}_n) P(y | \mathcal{D}_n)}{\gamma l(x | \mathcal{D}_n) + (1-\gamma)g(x|\mathcal{D}_n)}dy
     = \left( \gamma + \frac{g(x | \mathcal{D}_n)}{ l(x|\mathcal{D}_n) }(1-\gamma) \right)^{-1} \int_{-\infty}^{y^*} (y^* - y) P(y | \mathcal{D}_n) dy
     \propto \left( \gamma + \frac{g(x | \mathcal{D}_n)}{ l(x|\mathcal{D}_n) }(1-\gamma) \right)^{-1}


 \int_{-\infty}^{y^*} P(y|\mathcal{D}_n)dy = \gammaに注意。


以上の結果から、 \alpha_{EI}(x)を最大化するには \frac{g(x|\mathcal{D}_n)}{l(x|\mathcal{D}_n)}を小さくする点 xを考えれば良いと分かります。


特に l(x|\mathcal{D}_n)の定義から、以下の手続きで獲得関数を最大化する点が求まります。
1.  x \sim l(x|\mathcal{D}_n)を複数サンプリングする
2. サンプリングした xから、最も \frac{g(x | \mathcal{D}_n)}{l(x | \mathcal{D}_n)}を小さくする点 x \in \mathcal{X}を次の入力として採用

以上がEIとTPEを用いたベイズ最適化法の手続きとなります。


ガウス過程回帰と比べて計算量が小さく、観測した点数が多い場合にも使えるのがメリットです。
Optunaのdefaultの最適化アルゴリズムです。


こちらの記事が分かりやすいです。
gochikika.ntt.com





3-2. 獲得関数

獲得関数についてもまとめておきます。
統計モデルから求めた獲得関数を最大化する点 xを次の入力に用います。

EI(Expected Improvement)

EI(期待改善量)はblackbox最適化で用いられる最も代表的な獲得関数です。
単純かつ比較的高速に動作するので、手っ取り早く試すのに向いている手法です。

EIの獲得関数は、閾値 y^*に対して \alpha_{EI}(x) := \mathbb{E}[\max (0, y^* - y)]として定義されます。


 y^*には様々な決め方が存在しますが、 y_1,..., y_n \mathcal{D}_nに含まれる観測値として、
 y^* := \min \{y_1,...,y_n\}と定める方法が用いられたりします。



LCB(Lower Confidence Bound)

多碗bandit問題や強化学習でもおなじみの考え方です。

LCBは信頼区間の下限を用いた獲得関数であり、
GPRなど明示的に分散も求まる統計モデルを用いて以下の様に定義されます。

 \alpha_{LCB}(x) := - (m_n(x) - \sqrt{\beta_n} \sigma_n (x))
ただし、 \beta_nはハイパーパラメータで  \beta_n = c \log n\ \ (c:const.)などとして定めます。
(最大化問題を考える場合はUCBを考えれば良いです。)

LCBは"不確かな時は楽観的に"というヒューリスティックに基づいていて、
標準偏差を期待値に足し合わせた評価値で探索と活用のトレードオフを取ります。


直感的かつ単純な実装で実現できる事や、リグレット等の理論解析のしやすさがメリットとなりますが、
実用上性能がそこまで良くない事も実験的に知られています。



Thompson Sampling

Thompson Samplingは、目的関数の事後分布から関数を一つサンプリングし、
その関数を最小化する点を次の入力として採用する手法です。

こちらもLCB/UCBと合わせて多碗バンディット問題でよく使われます。
blackbox最適化は連続腕bandit問題とも解釈出来るので、banditの手法が幾つか登場します。


ガウス過程を用いる場合、Thompson Samplingの獲得関数は以下の様に定義されます。
 \alpha_{TS}(x) = - f^{(n)}(x),\ \ \ \ (f^{(n)} \sim \mathcal{GP}(m_n, k_n))


 f^{(n)}は確率的に定まる関数であり、直接 \alpha_{TS}(x)を最大化するのは難しいので、
今回は \mathcal{X}を有限個に分割して近似的に f^n(x)を計算します。


まず、 \mathcal{X} = \sum_{i=1}^\Delta I_iとして互いに素な集合族 \{I_i\}_{i=1}^\Deltaを取り、
各集合 I_iから代表点 p_i \in I_i \subset \mathcal{X}を取ります。
この代表点を用いて、 \mathcal{X} \simeq \{p_i\}_{i=1}^\Deltaとして有限個の点で \mathcal{X}を近似します。

さらに、 \mathcal{GP}(m_n, k_n)を用いて p_iに対応する q_i := f(p_i) \sim \mathcal{GP}(m_n, k_n)をサンプルし、
 \{(p_i, q_i) \}_{i=1}^\Deltaによって、関数 f^{(n)} f^{(n)} \simeq \hat{f}^{(n)}: p_i \mapsto q_iとして有限個の点で近似します。


これらの近似によって、以下の単純な計算で \alpha_{TS}(x)を最大化する点を近似的に求められる様になります。
 \mathrm{argmax}_{x \in \mathcal{X}}\ \alpha_{TS}(x) \simeq \mathrm{argmax}_{p_i \in \{x_1,...,x_\Delta\}} \left( - \hat{f}^{(n)}(p_i) \right) = \mathrm{argmin}_{i = \{1,...,\Delta\}}\ q_i

勿論この計算法では予め集合 \mathcal{X}の分割法を定める必要があります。


 \alpha_{TS}(x)を最大化する xは、 fを最小化する点の事後分布 p_{min}(x^* | \mathcal{D}_n)に従う点としても解釈できます。
この観点から捉えると、banditでのThompson Samplingとの関連性が分かりやすいと思います。



Thompson Samplingは多碗bandit問題では有力な手法なのですが、
ベイズ最適化ではそれほど使われていない様です。(高次元で性能が悪い)



※ banditの手法の比較は以下で試せます。(気になる方は試してみて下さい)
github.com




(補足)獲得関数の最小化について

特にGPRを用いた獲得関数の最小化について補足します。
先程定義した m_n(x) \sigma_n(x) = k_n(x, x) + \sigma^2を考えます。


これらの関数の x^{(i)}についての勾配は以下の様に求まります。( m_0(x) =0としています)
 \frac{\partial m_n (x)}{\partial x^{(i)}} = \left( \frac{\partial \mathbf{k}_n (x)}{\partial x^{(i)}} \right)^\top (\mathbf{K}_n + \sigma^2 I_n)^{-1} \mathbf{y}_n
 \frac{\partial \sigma_n(x)}{\partial x^{(i)}} = \frac{1}{2\sigma_n (x)} \left( \frac{\partial k_0(x, x)}{\partial x^{(i)}} - 2 \left( \frac{\partial \mathbf{k}_n(x)}{\partial x^{(i)}} \right)^\top (\mathbf{K}_n + \sigma^2 I_n)^{-1} \mathbf{k}_n(x)  \right)


これらの勾配は \mathbf{K}_n + \sigma^2 I_nを事前計算しておく事で、計算量を O(n^2)まで削減できます。

EIやLCBの獲得関数の勾配は m_n(x),\ \sigma_n(x)の勾配を用いて表現されるので、
上記の行列を事前計算しておく事で、獲得関数の最適化計算の計算量を削減できます。

今回の実装では事前計算を行っておりません






4. 実装(GP-LCB)

GPR+LCBを用いた場合のプログラムを実装してみます。
この手法はGP-LCB(GP-UCB)と呼ばれ、連続腕バンディット問題の方策の1つです。


ここからは、 f(x) = \exp(x/2) + \sin(2\pi x)を最小化する xblackbox最適化で求める事を考えます。
ただし \mathcal{X} = [-1.0, 1.0]とします。


GPRはscikit-learnのGaussianProcessRegressorを使用し、
kernel関数にはmatern kernelにwhite kernel、Constant kernelを組み合わせた物を用います。

matern kernelの定義は以下を参照して下さい。
sklearn.gaussian_process.kernels.Matern — scikit-learn 1.3.2 documentation


また、LCBの \beta_nには \beta_n := 2.0 * \log nを使用し、
 f(x)の観測ノイズ \epsilon \epsilon \sim \mathcal{N}(0, 0.01)とします。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import Matern, WhiteKernel, ConstantKernel
from scipy.optimize import minimize
from functools import partial


def f(x: np.ndarray) -> np.ndarray:
    x = x.reshape(-1)
    return np.exp(x/2.0) * np.sin(2 * np.pi * x)


def obs(value: np.ndarray, noise: bool | float = False) -> np.ndarray:
    return value + np.random.normal(0, noise, size=value.shape) if noise else value


def plot_data(xs: np.ndarray, ys: np.ndarray, x_bound: np.ndarray, gpr: GaussianProcessRegressor, n) -> None:
    x = np.linspace(x_bound[:, 0], x_bound[:, 1], 200).reshape(-1)
    y, sigma = gpr.predict(x.reshape(-1, 1), return_std=True)
    plt.plot(x, f(x), label='f(blackbox)')
    plt.plot(x, y, label='GPR(model)')
    plt.fill_between(x, y - sigma, y + sigma, alpha=0.2, label='sigma')
    plt.fill_between(x, y - sigma * np.sqrt(2.0 * np.log(n)), y + sigma * np.sqrt(2.0 * np.log(n)), alpha=0.2, label='UCB/LCB')
    plt.scatter(xs, ys, marker='x', color='red')
    plt.xlabel("x")
    plt.ylabel("f(x)")
    plt.legend()
    plt.show()


def lcb_policy(x_bound: np.ndarray, gpr: GaussianProcessRegressor, n: int) -> np.ndarray:
    def lcb_score(x: np.ndarray, gpr: GaussianProcessRegressor, n: int) -> float:
        mean, sigma = gpr.predict(x.reshape(1, -1), return_std=True)
        score = (mean - sigma * np.sqrt(2.0 * np.log(n)))[0]
        return score

    objective = partial(lcb_score, gpr=gpr, n=n)
    point_num = 30
    x_next_candidates: list = []

    for x in np.random.uniform(x_bound[:, 0], x_bound[:, 1], size=(point_num, x_bound.shape[0])):
        # 多点スタートで最適化 (local minimum対策)
        ans = minimize(objective, x0=x, bounds=x_bound)
        x_next_candidates.append([ans.fun, ans.x])

    x_next = min(x_next_candidates)[1]
    return x_next


def GP_LCB(x_bound: np.ndarray, trial_num: int):
    kernel = Matern(nu=1.5) + WhiteKernel() + ConstantKernel(1.0)  # カーネルの設定
    gpr = GaussianProcessRegressor(kernel=kernel, alpha=0, random_state=0)  # モデルの作成

    # observed data
    xs: np.ndarray = np.empty(0)
    ys: np.ndarray = np.empty(0)

    # 初期値の設定
    x_next: np.ndarray = np.random.uniform(x_bound[:, 0], x_bound[:, 1], size=(x_bound.shape[0],))

    for i in range(trial_num):
        print(f"trial: {i+1}")
        y_next = obs(f(x_next), noise=0.01)  # 観測値の取得
        xs = np.append(xs, x_next).reshape(-1, x_next.shape[0])
        ys = np.append(ys, y_next)
        gpr.fit(xs, ys)  # モデルの学習
        plot_data(xs, ys, x_bound, gpr, i+1)
        x_next = lcb_policy(x_bound, gpr, i + 1)

    return min(ys)


if __name__ == '__main__':
    x_bound_array = np.array([[-1.0, 1.0]])
    GP_LCB(x_bound_array, 20)

※ GPRのkernelの選び方などGP-LCB自体にもハイパラが多く、上手く最小化出来ない時もあります。


最適化計算等めんどうな部分はscikit-learnのモジュールで誤魔化しています。

本来GPRを用いた場合の獲得関数の最適化計算には逆行列の事前計算テクが使えるのですが、
今回は省略させて頂きました。


以下がプログラムの実行結果となります。
青の帯が分散を表し、オレンジの帯の上端がUCBを下端がLCBを表しています。

赤のxが入力した点 xを表しています。

反復に従って、 fを最小化する点が求まっている事が分かります。





5. optunaを用いた最適化

・optuna公式
optuna.org


5-1. install

以下のコマンドでinstall出来ます。

pip install optuna


dashboardで可視化したい人は以下もinstallして下さい。

pip install optuna-dashboard

installで躓く事は無いと思います。



5-2. toy-modelに対する最適化(sample code)

簡単な関数に対して、optunaを用いたblackbox最適化を実装してみます。


先ずは1変数で考えます。
最適化したい変数をtrial objectで記述し、以下のcodeで最適化します。

import optuna

def objective(trial):
    x = trial.suggest_float('x', -10, 10)  # xの範囲を[-10, 10]とする
    return (x - 2) ** 2

study = optuna.create_study()  # problem setting
study.optimize(objective, n_trials=100)  # blackbox最適化の実行
study.best_params  # 結果

結果はstudy objectのパラメータに格納されます。
optuna.study.Study — Optuna 3.5.0 documentation


また、optunaでは以下の様にアルゴリズムが選択出来ます。
Efficient Optimization Algorithms — Optuna 3.5.0 documentation

デフォルトではTPEを使用します。



2変数になっても同じ様に実装出来ます。

import optuna

def func(x, y):
    return x**2+y**2+x*y

def objective(trial):
    x = trial.suggest_float('x', -1, 1)  # xの範囲指定
    y = trial.suggest_float('y', -1, 1)  # yの範囲指定
    return func(x, y)

study = optuna.create_study()
study.optimize(objective, n_trials=100)
study.best_params

簡単にblackbox最適化を試す事が出来ます。



5-3. 結果の可視化

以下のコードで結果を可視化出来ます。

# Empirical Distribution Function Plot
fig = optuna.visualization.plot_edf(study)
fig.write_image('edf.png') 

# Optimization History Plot
fig = optuna.visualization.plot_optimization_history(study)
fig.write_image('optimization_history.png')

# parallel coordinate plot
fig = optuna.visualization.plot_parallel_coordinate(study)
fig.write_image('parallel_coordinate.png')

# slice
fig = optuna.visualization.plot_slice(study)
fig.write_image('slice.png') 


こちらの記事に詳しく説明されています。(コードもこちらの記事から引用しました。)
ibetoge.hatenablog.com




5-4. マルチプロセス化

optunaはRDBを用いたマルチプロセス化にも対応しています。
今回はSQLiteを用いて実装します。


以下のコマンドでデータベースを生成する事が出来ます。(pythonからでも生成出来ます)

optuna create-study --study 'optuna_test' --storage 'sqlite:///optuna_testdb.db'

pythonでも生成出来ます。

DATABASE_URI = 'sqlite:///optuna_testdb.db'
study_name = 'optuna_test'
optuna.create_study(study_name=study_name, storage=DATABASE_URI)

データベースが生成できたら、次はmulti-processでの実行を実装します。
今回はmultiprocessingを用いて実装します。


multi-processでoptunaを実行するコードは以下の様になります。
簡単かつ直感的に実装出来る様になっています。

import optuna
import multiprocessing
from multiprocessing import Process


def func(x, y):
    return x ** 2 + y ** 2 + x * y


def objective(trial):
    x = trial.suggest_float('x', -1, 1)  # xの範囲指定
    y = trial.suggest_float('y', -1, 1)  # yの範囲指定
    return func(x, y)


def optimize(study_name, storage, n_trials):
    study = optuna.create_study(study_name=study_name, storage=storage, load_if_exists=True)
    study.optimize(objective, n_trials=n_trials)


n_trials = 10
concurrency = multiprocessing.cpu_count()  # CPUの数を取得
n_trials_per_cpu = n_trials / concurrency

DATABASE_URI = 'sqlite:///optuna_testdb.db'
study_name = 'optuna_test'
optuna.create_study(study_name=study_name, storage=DATABASE_URI)  # dbの生成

# multi-process化
workers = [Process(target=optimize, args=(study_name, DATABASE_URI, n_trials_per_cpu)) for _ in range(concurrency)]

for worker in workers:
    worker.start()

for worker in workers:
    worker.join()

study = optuna.load_study(study_name=study_name, storage=DATABASE_URI)  # 結果の取得

fig = optuna.visualization.plot_edf(study)
fig.write_image('edf.png')

# Optimization History Plot
fig = optuna.visualization.plot_optimization_history(study)
fig.write_image('optimization_history.png')

# parallel coordinate plot
fig = optuna.visualization.plot_parallel_coordinate(study)
fig.write_image('parallel_coordinate.png')

# slice
fig = optuna.visualization.plot_slice(study)
fig.write_image('slice.png')


公式Documentにも記述があります。
Easy Parallelization — Optuna 3.5.0 documentation





6. 終わりに

ベイズ最適化以外の手法やアルゴリズムの使い分けにも触れようと思っていたのですが、
カレンダーの日にちまでに間に合いませんでした。

申し訳ないです...。






参考文献

この記事を書くに当たって参考にした解説/本を列挙しておきます。


blackbox最適化の幅広いサーベイ
search.ieice.org


ベイズ最適化
www.kindaikagaku.co.jp


ガウス過程回帰
www.jstage.jst.go.jp


・bandit問題
www.kspub.co.jp