【ガウス過程回帰】カーネルの紹介とボストン住宅価格予測

こんばんは。最近、ガウス過程回帰（Gaussian Process Regression, GPR）なる技術を勉強する機会がありました。誤差付きで予測できる便利な技術、と理解していますがGPRとニューラルネットワーク（NN）はどちらが予測精度が高いの？という素朴な疑問があったので予測精度の比較をやりました。参考書などを読むと、「単層のノード数無限ニューラルネットワークとGPRは等価である」などと書いてありますが、正直腑に落ちません。。。というわけで、理解を深める一助として様々なカーネルを使いながらGPRとNNを用いてボストン住宅価格データセットを用いて予測をしました。
〈目次〉

そもそもガウス過程回帰ってなあに？
ガウス過程回帰による予測
ニューラルネットワークによる予測
GPRの欠点と使いどころ

そもそもガウス過程回帰ってなあに？

ガウス過程回帰（GPR）の詳細な説明はここでは行いませんが、以下の特徴を備えた強力な回帰手法として知られています。

観測データに対して内挿を行う。
予測値だけでなく経験的信頼区間を同時に計算できる。つまり予測値の誤差を評価することができる。
問題によって適切なカーネルを定義することで、予測対象の滑らかさや周期性を上手く再現することができる。
似た入力データ同士は似た出力をする。

特にGPRを魅力的な手法足らしめているのは２と３だと私は考えています。予測誤差を計算できるため予測値がどの程度信頼できるかの判断材料になりますし、複数のカーネルを用いることで様々なデータ構造に柔軟に対応できます。
なお、ここまで気軽にカーネルカーネル言っていますが、その理論的深みを私は理解していません。ここでは深入りせず、入力データ点ごとの共分散であると大雑把にお考え下さい。詳しくは以下の本や、Scikit-learnのウェブページをご参照ください。

ガウス過程と機械学習 (機械学習プロフェッショナルシリーズ)

作者:持橋大地,大羽成征
発売日: 2019/03/09
メディア: 単行本（ソフトカバー）

scikit-learn.org

例として、ガウス過程回帰を使うと以下のような図が得られます。これはsin関数をガウス過程回帰で内挿したものですが、予測値に加えて、標準偏差や95%信頼区間を合わせて表示させることができます。

f:id:Dajiro:20200421002619p:plain — ガウス過程回帰。真値、予測値、95%信頼区間を表示。

ガウス過程回帰による予測

それでは、ガウス過程回帰をさっそく使ってみましょう。実装はscikit learnのモジュールを使って簡単にできます。scikit learnに備え付けのカーネルや自作のカーネルを自由に使うことができ、ハイパーパラメータは初期値さえ設定すれば勾配上昇法により自動的に計算されます（便利ですねぇ）。逆に細かいことを知らなくても結果が出てしまうという点には注意しなければいけません。
今回の計算では以下の３つのカーネルを用いて計算しました。

Radial-basis function (RBF)カーネル： $\rm{exp}$ $(-\frac{1}{2}d(\frac{xi}{l}, \frac{xj}{l})^{2})$

RBFカーネルは長さスケールを調整するハイパーパラメータが含まれていますが、これは最尤推定により最適化されます。またRBFカーネルを用いることで予測対象の滑らかさを表現することができます。

Rational quadratic カーネル： $(1 + \frac{d(xi, xj)^{2}}{2\alpha l^{2}})^{-\alpha}$

Rational quadratic カーネルは異なる長さスケールを持つRBFカーネルの無限和としてとらえられ、α→∞とするとRBFカーネルと一致する性質を持つそうです。２つのハイパーパラメータが最適化されます。

Matern カーネル：(長いので数式省略・・・以下リンクをご参照ください！)

1.7. Gaussian Processes — scikit-learn 0.23.1 documentation
RBFカーネルの一般化と言えるカーネルです。νという新しいパラメータが付け加わっており、νが∞ならRBFカーネルと一致し、v = 1/2, 3/2, 5/2の場合それぞれ異なる性質をもったカーネルに変化します。これは自動決定されず、ユーザーが決めるパラメータになります。

本来は予測対象のデータ構造を分析した上で適切なカーネルを設計すべきなのでしょうが、とりあえず代表的な上記３つのカーネルを使って計算してみることにしましょう。コードは以下の通りです。さぁ、回帰してみましょう。

#もろもろの便利ツールをインポート
import numpy as np
import seaborn
from sklearn.model_selection import train_test_split
from sklearn.datasets import load_boston
import pandas as pd
import matplotlib.pyplot as plt
plt.style.use('seaborn-muted')
import japanize_matplotlib
#ガウス過程回帰用のモジュール
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import ConstantKernel, RBF, WhiteKernel, ExpSineSquared, Matern, RationalQuadratic
from sklearn.preprocessing import StandardScaler

# 住宅価格の読み込みとtrain, testの分割
dataset = load_boston()
y = dataset.target
X = dataset.data
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.20,  random_state=42)
_, n_feat = X_train.shape

#標準化
X_mean = {i: np.mean(X_train[:, i]) for i in range(n_feat)}
X_std  = {i: np.std(X_train[:, i]) for i in range(n_feat)}
y_mean = np.mean(y_train)
y_std  = np.std(y_train)

for i in range(n_feat):
    X_train[:, i] = (X_train[:, i] - X_mean[i]) / X_std[i]
    X_test[:, i]  = (X_test[:, i] - X_mean[i]) / X_std[i]
y_train[:] = (y_train[:] - y_mean) / y_std

# GPモデルの構築
kernels = [1.0 * RBF() + WhiteKernel(),
          1.0 * RationalQuadratic(length_scale=1.0, alpha=0.1) + WhiteKernel(),
          1.0 * Matern(length_scale=1.0, length_scale_bounds=(1e-1, 10.0),
                        nu=1.5) + WhiteKernel()]
#３つのカーネルに対してループを回し回帰と結果のプロットを実行
for n, kernel in enumerate(kernels):
    gpr = GaussianProcessRegressor(kernel=kernel, alpha=0)
    gpr.fit(X_train, y_train)

    print(gpr.kernel_)
    pred_mu, pred_sigma = gpr.predict(X_test, return_std=True)
    pred_mu = pred_mu * y_std + y_mean
    pred_sigma = pred_sigma.reshape(-1, 1) * y_std
    print(np.corrcoef(pred_mu, y_test)[0][1])
    fig = plt.figure(figsize=(6, 6))
    if n == 0:plt.title('RBF kernel', fontsize = 16)
    if n == 1:plt.title('Rational quadratic kernel', fontsize = 16)
    if n == 2:plt.title('Matern kernel', fontsize = 16)
    plt.errorbar(y_test, pred_mu, yerr = pred_sigma, fmt = "o")
    plt.tick_params(labelsize=16)
    x = np.linspace(0, 50)
    y = np.linspace(0, 50)
    plt.plot(x, y, 'k--')
    plt.xlabel("真値", fontsize = 16)
    plt.ylabel("予測値", fontsize = 16)
    plt.grid(alpha = 0.5)
    plt.savefig(f"GP_boston_{n}.png")
    plt.show()

すると以下のような図が得られます。相関係数はRBF, Rational quadratic, Maternカーネルでそれぞれ相関係数は0.921, 0.924, 0.925程度の予測精度です。精度自体は、この問題に関しては大きな差はないという結果ですね。特筆すべきは標準偏差が得られるのでエラーバーを付けて可視化できることです。より詳細な予測精度を見積もりたい場合にはNNよりこちらの方が望ましいのではないでしょうか。

f:id:Dajiro:20200428094906p:plain — RBF, RationalQuadratic, Maternカーネルの順に結果をプロット（誤差は標準偏差σ）。

f:id:Dajiro:20200428094916p:plain — RBF, RationalQuadratic, Maternカーネルの順に結果をプロット（誤差は標準偏差σ）。

ニューラルネットワークによる予測

続いて、ニューラルネットワークで予測してみましょう。同じくボストンの住宅価格予測を行います。以下のようなコードを書きました。本当はハイパーパラメータ最適化すべきなのですが、ここでは私が手動で設定したものを使っています。隠れ層数は3で、ノード数は100、活性化関数はreluとしました。

#ニューラルネットワーク回帰用のkerasモジュールをインポート
from tensorflow.python import keras as K

# 住宅価格の読み込みとtrain, testの分割と標準化はガウス過程回帰と同じ
（省略）

#モデルの定義
model = K.Sequential([
    K.layers.Dense(100, input_shape=(n_feat,), activation="relu"),
    K.layers.Dense(units=100, activation="relu"),
    K.layers.Dense(units=100, activation="relu"),
    K.layers.Dense(units=1)
])

#モデルのコンパイルと学習実行
model.compile(loss="mean_squared_error", optimizer="Adam")
model.fit(X_train, y_train, epochs=100)

# 予測
for i in range(n_feat):
    X_test[:, i] = (X_test[:, i] - X_mean[i]) / X_std[i]
predicts = model.predict(X_test)

f:id:Dajiro:20200428095621p:plain — ニューラルネットワークによる予測値。

相関係数を計算すると0.930という精度が得られました。単に精度を比較するならNNがGPRより高いわけですが、それもわずかな差ですしハイパーパラメータも最適化していないので比較にはあまり意味がないでしょう。少なくとも今回の計算条件では同等の精度が得られることが可視化されました。

GPRの欠点と使いどころ

GPRとNNはどのように使い分けすべきでしょうか。GPRは予測誤差を出力できる素晴らしい方法ですが、以下のような欠点があります。

データ数が増大すると計算効率が非常に悪くなる。

重要な入力のみに着目するということはなく（スパースではなく）、入力全てを活用する。

画像やグラフ構造のデータ処理には向かない。また出力も基本的には1チャネル。

これらの制約がかからない問題設定ではGPRは積極的に活用すべきだと思います。また確率的な分類問題にも適用できるので、無茶なことを言わなければ非常に有用な飛び道具と言えるでしょう。逆に、画像、音声、言語、グラフなど様々な問題に適用でき、人間の無茶な要求にもこたえ続けているニューラルネットワークが異常な存在とも言えますが・・・

おしまい。