ハイパーパラメータチューニングにあたっては，一般的にグリッドサーチなどの方法が用いられます．しかし，一回の学習に時間がかかるようなモデルにおいては，より効率的なハイパーパラメータ探索が求められます．
今回はガウス過程回帰を利用することによって，より効率的なハイパーパラメータチューニングを目指します．

ガウス過程回帰によるハイパーパラメータ調整

獲得関数をUCB(Upper Confidence Bound)とし，ガウス過程回帰に基づいたハイパーパラメータチューニングの方法をまとめます．
アルゴリズムの構築にあたっては，以下の大変面白い記事を参照しています．

問題設定

訓練データ$(X_{\rm train}, {\bf y}_{\rm train})$を用いてモデルを作成し，テストデータの入力$X _ {\rm test}$に対して出力${\bf y} _ {\rm test}$を予測する問題を考えます．今回のモデルでは，ハイパーパラメータ${\bf \theta}$のもとで学習したパラメータ${\bf w} ^ { * }$に基づいて，予測値${\bf \hat{y}} _ {\rm test}$を出力するものとします．ここで，予測の精度を測定するにあたって何かしらのスコア$s = s({\bf \hat{y}} _ {\rm test}, {\bf y} _ {\rm test})$を用います．本記事でのハイパーパラメータチューニングは，スコア$s$を最大化するようなハイパーパラメータ${\bf \theta}^*$を求めることであると定義します．

最適化対象の獲得関数

スコア$s$を最大化することが目的ではあるものの，ハイパーパラメータ${\bf \theta}$の空間$\Theta$全体の走査は(特に次元が大きい場合)難しいです．そこで，①学習の初期段階においては探索が進んでいない領域を積極的に探し，②学習が進んだ段階では高いスコアが得られた近辺を捜索する，というアルゴリズムを考えるため，以下のような定式化を考えます．

まず，$n$個のハイパーパラメータの組$({\bf \theta}_1, \cdots, {\bf \theta}_n)$に対して，スコアが${\bf s} = (s_1, \cdots, s_n)$のように得られているものとします．ガウス過程回帰を用いると，新しいハイパーパラメータ${\bf \theta} _ {n+1}$に対するスコアの予測値$\hat{s} ({\bf \theta} _ {n+1})$とその標準偏差$\sigma({\bf \theta} _ {n+1})$が求められます．
獲得関数をUCB(Upper Confidence Bound)とすると，以下を最大化する${\bf \theta} _ {n+1}$を探索することになります：

\alpha({\bf \theta} _ {n+1}) = \hat{s} ({\bf \theta} _ {n+1}) + \sqrt{\dfrac{\log (n+1)}{n+1} } \sigma({\bf \theta} _ {n+1})

第一項はスコアが高いところを，第二項は標準偏差が大きいところ(未探索の箇所)を探索するインセンティブとなります．特に後者の効果は，$n$が小さいほど(探索が進んでいない学習の初期において)大きくなります．

アルゴリズム

実際にはハイパーパラメータ全体の空間$\Theta$全体を対象とできないため，適当な組を指して$\Theta$で表記することにします．またガウス過程回帰のため，あらかじめ抽出した$n_0$個のハイパーパラメータ$\Theta _ {n _ 0} = ({\bf \theta}_1, \cdots, {\bf \theta} _ { n _ 0 } )$に対して，スコアが${\bf s} = (s_1, \cdots, s _ { n _ 0 } )$として求められているものとします．
今回のアルゴリズムは，以下のような手順で与えられます：

$n$個のハイパーパラメータ$\Theta _ {n}$を入力，対応するスコア${\bf s}$を出力として，ガウス過程回帰．
ハイパーパラメータの空間$\Theta$の中で未探索の組について獲得関数$\alpha({\bf \theta} _ {n+1})$を計算し，UCBを最大化させる${\bf \theta} _ {n+1}$を求める．
求めた${\bf \theta} _ {n+1}$をもとに作成したモデルをもとに${\bf y} _ {\rm test}$の予測値${\bf \hat{y}} _ {\rm test}$に対してスコア$s_{n+1}$を算出．

数値実験

適当に作成したデータに対して，ベイズ線形回帰による予測を目指します．ベイズ線形回帰を用いたのは，後述のようにエビデンス$p({\bf t} \mid \alpha, \beta)$を明示的に計算できることが理由です．ベイズ線形回帰の実装にあたっては，以下の記事を参考にしています．

回帰モデル

$N$個の訓練データ$x_1, \cdots, x_N$に対し，出力${\bf t} = (t_1, \cdots, t_N)^t$が得られているものとします．
いま$n$番目の訓練データ$x_n$を$d$次元の特徴量ベクトル${\bf \phi}_n = (1, x_n^1, \cdots, x_n^{d-1})^t$に変換し，以下の特徴量行列$\Phi$を作成します．

\begin{eqnarray}
\Phi = \left(
\begin{array}{c}
{\bf \phi}_1^t \\
\vdots \\
{\bf \phi}_N^t
\end{array}
\right)
\end{eqnarray}

以下のベイズ線形回帰モデルを仮定します．

p({\bf t} \mid {\bf w}, \beta) = \left( \dfrac{\beta}{2 \pi} \right)^{N/2} \exp \left[ -\dfrac{\beta}{2} ({\bf t} - {\bf w}^t \Phi)^2 \right] \\
p({\bf w} \mid \alpha) = \left( \dfrac{\alpha}{2 \pi} \right)^{d/2} \exp \left[ -\dfrac{\alpha}{2} {\bf w}^2 \right]

ただし${\bf w}$をパラメータとし，$\alpha, \beta$をハイパーパラメータとしています．このとき，${\bf w}$の事後分布$p({\bf w} \mid {\bf t})$は

p({\bf w} \mid {\bf t}) = \mathcal{N}({\bf w} \mid {\bf m}_N, S_N) \\
{\bf m}_N = \beta S_N \Phi^t {\bf t} \\
S_N = (\alpha 1_d + \beta \Phi^t \Phi)

として計算されます．
${\bf w}$を積分消去することにより，ベイズ線形回帰の場合はエビデンス$p({\bf t} \mid \alpha, \beta)$を明示的に計算できます：

\log p({\bf t} \mid \alpha, \beta) = \dfrac{d}{2} \log \alpha + \dfrac{N}{2} \log \beta - \dfrac{\beta}{2} |{\bf t} - \Phi {\bf m}_N |^2 - \dfrac{\alpha}{2} {\bf m}_N^2 - \dfrac{1}{2} \log | \alpha 1_d + \beta \Phi^t \Phi | - \dfrac{N}{2} \log 2 \pi

具体的なモデル

最適化したいモデルとしては，以下を想定します：

t = \sin(x) + 0.5 \sin(2x) + 0.1x + 0.4 \epsilon \\
\epsilon \sim \mathcal{N}(0, 1)

$x$の値域は$x \in [-\pi, \pi]$であり，一様分布からサンプルするものとします．

$N = 50$個のデータをプロットしたものが下図です．

明示的に計算したエビデンス$p({\bf t} \mid \alpha, \beta)$を等高線でプロットしたものが左図，エビデンス最大化によって最適化したモデルが右図です．予測したモデルは，ある程度真のモデルと近い形になっていることが確認できます．

ハイパーパラメータチューニングの実装

一般的にはエビデンス$p({\bf t} \mid \alpha, \beta)$を厳密に計算することはできないため，スコアとして決定係数を利用します．いくつか定義されていない変数がありますが，おおよそは以下のように実装できます．全体の実装については，私のgithubリンクから確認できます．

import itertools # ハイパーパラメータの組み合わせを生成するために使用
from sklearn.gaussian_process import *

# ガウス過程回帰の実装
kernel = kernels.RBF(1.0, (3e-1, 1e3)) + kernels.ConstantKernel(1.0, (3e-1, 1e3)) + kernels.WhiteKernel()
model_hyper_param = GaussianProcessRegressor(
    kernel = kernel,
    alpha = 1e-3,
    optimizer = "fmin_l_bfgs_b",
    n_restarts_optimizer = 20,
    normalize_y = True)

# スコアを算出しているdf_params_w_scoreをもとにガウス過程回帰し，df_to_append_ucbにUCBを追加する(対数変換に注意)
def upper_confidence_bound(df_params_w_score, df_to_append_ucb):
    # Gauss過程回帰の学習
    x_hyper_param_train = np.log(df_params_w_score.values[:, :-1]) # スコア以外を対数変換
    y_hyper_param_train = df_params_w_score[["score"]].values
    model_hyper_param.fit(x_hyper_param_train, y_hyper_param_train) # ガウス過程回帰
    # Gauss過程回帰の予測
    x_hyper_param_test = np.log(df_to_append_ucb.values[:, :-1]) # スコア以外を対数変換
    score_pred, score_std = model_hyper_param.predict(x_hyper_param_test, return_std = True)
    score_pred = score_pred[:, 0] # 変な形になっているため，修正
    # スコアがわかっていないサンプルに対してUCBをくっつける
    df_to_append_ucb["UCB"] = score_pred + np.sqrt(np.log(len(df_params_w_score))/len(df_params_w_score)) * score_std
    return df_to_append_ucb

# ハイパーパラメータのチューニング
n_sample = 50 # 探索するハイパーパラメータの数
n_sample_for_GR = 3 # Gaussian Regressionのために必要とする，事前のハイパーパラメータ探索数

df_params = pd.DataFrame(list(itertools.product(alpha_list, beta_list)), columns = ["alpha", "beta"])
df_params["score"] = np.nan

for s in range(n_sample + 1):
    print("sample = {}".format(s), end = "\r")
    idx_score_nan = df_params[df_params["score"].isnull()].index # 未探索の箇所
    
    if s < n_sample_for_GR: # 試行回数がGauss過程回帰に十分でないとき
        idx_sample = np.random.choice(idx_score_nan)
    else:
        # スコアが分かっているサンプルとわかっていないものに分割
        df_params_w_score = df_params[~(df_params["score"].isnull())].copy() # スコアがわかっているサンプル
        df_params_wo_score = df_params[df_params["score"].isnull()].copy() # スコアがわかっていないサンプル
        df_params_wo_score = upper_confidence_bound(df_params_w_score, df_params_wo_score) # UCBを追加
        idx_sample = df_params_wo_score.sort_values(by = "UCB", ascending = False).index.values[0]
    
    # 探索すべきパラメータ
    alpha, beta, _ = df_params.iloc[idx_sample].values # scoreはいらないので捨てる
    model = Bayesian_Linear_Regression(alpha = alpha, beta = beta) # ベイズ線形回帰
    model.calc_posterior_params(X_train, t_train) # ハイパーパラメータをもとに最適化
    df_params.at[idx_sample, "score"] = model.calc_score(X_test, t_test) # 予測値に基づいてスコア算出

以下は，ガウス過程回帰によりハイパーパラメータを探索している際の獲得関数のプロットです．初期段階では全体を広く調べているものの，次第にスコアが高いところにあたりを付けて集中的に調べていることが確認できます．

まとめ

今回用いたベイズ線形回帰であれば，学習に対して時間はかからず，またエビデンスを明示的に計算することも可能です．一方で，学習に時間がかかるようなモデルであれば，ハイパーパラメータチューニングの効率化は極めて効果的です．特にチューニング部分は一般的な形で記述しているため，広いモデルに対しても利用可能であると考えています．
一方で，ガウス過程回帰を用いている都合上，ハイパーパラメータの空間が離散的であるような場合には取り扱いは難しくなります．あくまで連続的なハイパーパラメータに限って，チューニングが可能になるという理解が良いかと思われます．