📈

[NumPy] 簡単なガウス過程回帰モデルをゼロから実装する (ハイパーパラメータチューニングなしver.)

2022/08/12に公開約15,700字1件のコメント

次回: https://zenn.dev/mory22k/articles/e3222d02724251
前回: https://zenn.dev/mory22k/articles/429ce4ad3208d7

簡易的なガウス過程回帰モデルをnumpyで実装します。

import numpy as np
import matplotlib.pyplot as plt

元の関数

元の関数を設定

元の関数として、適当な1変数関数を設定しておきます。

f(x) = 3\sin(x) + 4\cos(2x) + \frac{1}{9}\exp\!\left(\frac{1}{3}x\right)
x_min = 0
x_max = 10

def func(x):
    return 3 * np.sin(x) + 4 * np.cos(2 * x) + 1/9 * np.exp(1/3 * x)

x = np.linspace(x_min, x_max)
y = func(x)

plt.figure(figsize=(12, 6))
plt.plot(x, y, color='black', linestyle='dashed', label='true')
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.legend()
plt.show()

output

データを得る

設定した元の関数からデータを得て、これを訓練データX_\textrm{train}として利用します。とりあえず、等間隔に15点クエリします。

num_init = 15

x_train = np.linspace(x_min, x_max, num_init)
y_train = func(x_train)

plt.figure(figsize=(12, 6))
plt.plot(x, y, label='true', linestyle='dashed', color='black')
plt.scatter(x_train, y_train, label='data', color='black', marker='x', s=200)
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.legend()
plt.show()

output

ブラックボックス関数とみなす

以降、設計した関数f(x)をブラックボックス関数として、その正体がわからないものとします。つまり、我々は次のようなものを見ているとイメージします。

output

ガウス過程回帰モデル

ガウス過程回帰モデルは、ベイズ線形回帰モデルをより一般化した回帰モデルです。ベイズ線形回帰モデルからの具体的な導出については、以下の記事をお読みください。

https://zenn.dev/mory22k/articles/2879abc4ae2db9

ガウス過程について

回帰関数\hat fは、平均関数が0、カーネル関数がk(x, x')のガウス過程に従うものと考えます。すなわち、

\hat f \sim \mathcal{GP} (0, k(x,x'))

です。このような場合、任意のn個の点\{x_1, x_2, \cdots, x_n\}に対して、出力\{f(x_1), f(x_2), \cdots, f(x_n)\}を積み重ねたベクトル

\bm f_n = \begin{bmatrix} f(x_1) \\ f(x_2) \\ \vdots \\ f(x_n) \end{bmatrix}

は、多変量ガウス分布に従います。

\bm f_n \sim \mathcal N(\bm 0, \bm K)
\bm K = \begin{bmatrix} k(x_1, x_1) & k(x_1, x_2) & \cdots & k(x_1, x_n) \\ k(x_2, x_1) & k(x_2, x_2) & \cdots & k(x_2, x_n) \\ \vdots & \vdots & \ddots & \vdots \\ k(x_n, x_1) & k(x_n, x_2) & \cdots & k(x_n, x_n) \end{bmatrix}

見方を変えて、\{x_1, x_2, \cdots, x_n\}\{f(x_1), f(x_2), \cdots, f(x_n)\}を訓練データと見なします。そのうえで、新たな点x_\astを入力した際の出力f_\ast := f(x_\ast)は、次のように予測することができます。

f_\ast | \bm f_n \sim \mathcal N\!\left( m, v \right)
\begin{align*} m &= \bm k_\ast^\mathsf{T} \bm K^{-1} \bm f_n \\ v &= k_{\ast\ast} - \bm k_\ast^\mathsf{T} \bm K^{-1} \bm k_\ast \end{align*}

ここで、

\begin{align*} \bm k_\ast &= \begin{bmatrix} k(x_1, x_\ast) \\ k(x_2, x_\ast) \\ \vdots \\ k(x_n, x_\ast) \end{bmatrix} \\ k_{\ast\ast} &= k(x_\ast, x_\ast) \end{align*}

です。

以上のものを実装する際は、具体的には次の2つを行います。

  1. カーネル関数k(x, x')の選択・設計
  2. \bm K, \bm k_\ast, k_{\ast\ast}の実装

カーネル関数の選択

カーネル関数k(x_p, x_q)は、共分散行列\bm K(p,q)成分に相当します。共分散行列の(p,q)行列は、x_px_qの共分散\textrm{Cov}[x_p, x_q]を表します。すなわち、

k(x_p, x_q) = \textrm{Cov}[x_p, x_q]

です。共分散\textrm{Cov}[x_p,x_q]はふつう、2つの確率変数x_p, x_qが「似ている」ほど値が大きくなります。したがって、通常、k(x_p,x_q)2つの入力x_p,x_qが近い位置にあるほど値が大きくなるように設計します。

世界中で、様々なカーネル関数が設計されています。ここでは、最もよく使われている動径基底関数カーネル (RBF Kernel; 二乗指数カーネル SE Kernel; あるいはガウスカーネル Gaussian Kernel)と呼ばれる、次のカーネル関数を実装します。

k_\textrm{RBF}(x, x') = \exp\!\left( -\frac{1}{\theta_0}\|x - x'\|^2 \right) \\ \theta_0 \gt 0
class rbf:
    def __init__(self, params):
        self.params = params
    
    def __call__(self, xi=0, xj=0):
        return np.exp(-1 / self.params[0] * (xi - xj)**2)

x'0に固定して、k_\textrm{RBF}(x, 0)をプロットしてみましょう。

kernel = rbf([1.])

x = np.linspace(-3, 3, 1000)
y = kernel(x, 0)

plt.figure(figsize=(12, 6))
plt.plot(x, y, color='#990033', label='RBF Kernel')
plt.xlabel('x')
plt.ylabel('k(x)')
plt.grid()
plt.legend()
plt.show()

output

x0を平均値とするガウス分布のような形状を描いていることがわかると思います。

続いて、パラメータの値を\theta_0 = 0.1に設定して、k_\textrm{RBF}(x, 0)を描いてみます。

kernel.params[0] = 0.1

x = np.linspace(-3, 3, 1000)
y = kernel(x, 0)

plt.figure(figsize=(12, 6))
plt.plot(x, y, color='#990033', label='RBF Kernel')
plt.xlabel('x')
plt.ylabel('k(x)')
plt.grid()
plt.legend()
plt.show()

output

\theta_0 = 1.0の場合に比べて、関数が尖った形状をしています。このように、パラメータを変更することでカーネルの形状を変更させることができます。これは、より正確な回帰分析を行う上でとても重要な性質になります。

カーネル行列の実装

続いて、\bm K, \bm k_\ast, k_{\ast\ast}を計算するための関数を実装します。

def kernel_matrix(xi, xj, kernel):
    xis, xjs = np.meshgrid(xj, xi)
    return kernel(xis, xjs)

ここで、訓練データの点数がnであるとき、各行列は、

\begin{align*} \bm K &\in \mathbb R^{n \times n} \\ \bm k_\ast &\in \mathbb R^{n \times 1} \\ k_{\ast\ast} &\in \mathbb R^{1 \times 1} \end{align*}

でなくてはならないことに注意しましょう。 点x_\astを適当に定めて計算し、形状を確認してみます。

kernel = rbf([1.])

K = kernel_matrix(x_train, x_train, kernel)

x_ = 1
k_ = kernel_matrix(x_train, x_, kernel)
k__ = kernel_matrix(x_, x_, kernel)

print(K.shape, k_.shape, k__.shape)
# (15, 15) (15, 1) (1, 1)

きちんと行列の形状の条件を満たしていることが確認できましたね。

予測分布を計算する

ある1つの点x_*に対する出力y_\ast := f(x_\ast)の予測分布の平均値mと分散vは、

\begin{align*} m &= \bm k_\ast^\mathsf{T} \bm K^{-1} \bm y_n \\ v &= k_{\ast\ast} - \bm k_\ast^\mathsf{T} \bm K^{-1} \bm k_\ast \end{align*}

によって計算することができます。試しに、x_\ast = 1として、y_\astの予測分布の平均値と分散を計算してみましょう。

kernel = rbf([1.])

K = kernel_matrix(x_train, x_train, kernel)
K_inv = np.linalg.inv(K)

x_ = 1
k_ = kernel_matrix(x_train, x_, kernel)
k__ = kernel_matrix(x_, x_, kernel)

k_K = k_.T.dot(K_inv)
y_mean = k_K.dot(y_train)
y_var = k__ - k_K.dot(k_)

print(y_mean, y_var)
# [1.08995174] [[0.00740021]]

同様の計算を、ある範囲x_\textrm{min} \le x \le x_\textrm{max}の点ひとつひとつに対して行い、グラフにプロットしてみましょう。これにより、その範囲の予測分布を図示することができます。

def predict(xs, x_train, y_train, kernel):
    K = kernel_matrix(x_train, x_train, kernel)
    K_inv = np.linalg.inv(K)

    y_mean = np.array([])
    y_var = np.array([])

    for x in xs:
        k_ = kernel_matrix(x_train, x, kernel)
        k__ = kernel_matrix(x, x, kernel)

        k_K = k_.T.dot(K_inv)
        y_mean = np.append(y_mean, k_K.dot(y_train))
        y_var = np.append(y_var, k__ - k_K.dot(k_))
    return y_mean, y_var

# 真の関数値
x = np.linspace(x_min, x_max)
y = func(x)

# 訓練データ
x_train = np.linspace(x_min, x_max, num_init)
y_train = func(x_train)

# 予測
kernel = rbf([1.])
y_mean, y_var = predict(x, x_train, y_train, kernel)
y_std = np.sqrt(np.maximum(y_var, 0))

# プロット
plt.figure(figsize=(12, 6))
plt.plot(x, y, label='true', linestyle='dashed', color='black')
plt.plot(x, y_mean, label='mean', color='#0066cc')
plt.fill_between(x, y_mean-y_std, y_mean+y_std, label='2σ credible region', color='#0066cc', alpha=0.3)
plt.scatter(x_train, y_train, label='data', color='black', marker='x', s=200)
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.legend()
plt.show()

output

元の関数をかなり上手く近似できていることがわかると思います。

予測分布計算時の注意

カーネル行列の正則性

さて、こうしてガウス過程回帰が行えるようになったわけですが、1つ注意点があります。次の例を見てみましょう。

x_train = np.array([0, 0])
y_train = func(x_train)

kernel = rbf([1.])
predict([0], x_train, y_train, kernel)
# LinAlgError: Singular matrix

訓練データX_\textrm{train}に全く同じ点が2つ以上存在すると、カーネル行列\bm Kの中に、全く同じ要素を持つ行が2つ生じてしまいます。例えば、X_\textrm{train} = \{x_1, x_2\}で、x_1 = x_2のような場合、

\bm K = \begin{bmatrix} k(x_1, x_1) & k(x_1, x_2) \\ k(x_2, x_1) & k(x_2, x_2) \end{bmatrix} = \begin{bmatrix} k(x_1, x_1) & k(x_1, x_1) \\ k(x_1, x_1) & k(x_1, x_1) \end{bmatrix}

であり、1行目と2行目が完全に一致してしまいます。それゆえ、\bm Kが正則行列でなくなってしまう (すなわち、特異行列: singular matrix になってしまう) ため、逆行列\bm K^{-1}が計算できなくなり、Singular matrixエラーが生じます。

ノイズを仮定する

このような事態を防ぐ簡単な方法の1つに、関数fを評価する際にガウスノイズが生じると仮定する方法があります。

y = f(x) + \varepsilon \\ \varepsilon \sim \mathcal N(0, \sigma^2)

これをまとめて書くと

y \sim \mathcal N(f(x), \sigma^2)

です。

ノイズを含むガウス過程

ノイズなしの場合には、任意のn個の点\{x_1, x_2, \cdots, x_n\}に対して、出力を積み重ねたベクトル

\bm f_n = \begin{bmatrix} f(x_1) \\ f(x_2) \\ \vdots \\ f(x_n) \end{bmatrix}

は、次のようなガウス分布に従うのでした。

\bm f_n \sim \mathcal N(\bm 0, \bm K)
\bm K = \begin{bmatrix} k(x_1, x_1) & k(x_1, x_2) & \cdots & k(x_1, x_n) \\ k(x_2, x_1) & k(x_2, x_2) & \cdots & k(x_2, x_n) \\ \vdots & \vdots & \ddots & \vdots \\ k(x_n, x_1) & k(x_n, x_2) & \cdots & k(x_n, x_n) \end{bmatrix}

一方、出力にノイズが生じることを仮定する場合には、

\bm y_n = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix} = \begin{bmatrix} f(x_1) + \varepsilon_1 \\ f(x_2) + \varepsilon_2 \\ \vdots \\ f(x_n) + \varepsilon_n \end{bmatrix} = \bm f_n + \bm \varepsilon_n

です。これは、多変量ガウス分布として

\bm y_n \sim \mathcal N(\bm f_n, \sigma^2 \bm I)

と書くことができます。これと\bm f_n \sim \mathcal N(\bm 0, \bm K)を合わせると、次のように書けます。

\bm y_n \sim \mathcal N(\bm 0, \bm K + \sigma^2 \bm I)

ノイズを含む予測分布

関数値の予測

\{x_1, x_2, \cdots, x_n\}\{y_1, y_2, \cdots, y_n\}を訓練データと見なします。そのうえで、新たな点x_\astを入力した際の出力f_\ast := f(x_\ast)は、次のように予測することができます。

f_\ast | \bm y_n \sim \mathcal N\!\left( m, v \right)
\begin{align*} m &= \bm k_\ast^\mathsf{T} (\bm K + \sigma^2 \bm I)^{-1} \bm y_n \\ v &= k_{\ast\ast} - \bm k_\ast^\mathsf{T} (\bm K + \sigma^2 \bm I)^{-1} \bm k_\ast \end{align*}

観測値の予測

観測値y_\astは、

y_\ast = f(x_\ast) + \varepsilon

から、次のように予測できます。

y_\ast | \bm y_n \sim \mathcal N\!\left( m, v_y \right)
\begin{align*} m &= \bm k_\ast^\mathsf{T} (\bm K + \sigma^2 \bm I)^{-1} \bm y_n \\ v_y &= k_{\ast\ast} - \bm k_\ast^\mathsf{T} (\bm K + \sigma^2 \bm I)^{-1} \bm k_\ast + \sigma^2 \end{align*}

ノイズを含むガウス過程回帰モデル

以上の考察をもとに、predict()関数を修正してみましょう。

def predict2(xs, x_train, y_train, kernel, sigma=1e-2, output_noise=True):
    K = kernel_matrix(x_train, x_train, kernel) \
        + sigma**2 * np.identity(len(x_train))
    K_inv = np.linalg.inv(K)

    y_mean = np.array([])
    y_var = np.array([])

    for x in xs:
        k_ = kernel_matrix(x_train, x, kernel)
        k__ = kernel_matrix(x, x, kernel)

        k_K = k_.T.dot(K_inv)
        y_mean = np.append(y_mean, k_K.dot(y_train))
        y_var = np.append(y_var, k__ - k_K.dot(k_) + sigma**2 * output_noise)
    return y_mean, y_var

output_noiseTrueの場合に観測値の予測を行い、Falseの場合に関数値の予測を行います。これを使って先ほどと同じ操作を行ってみます。

x_train = np.array([0, 0])
y_train = func(x_train)

kernel = rbf([1.])
y_mean, y_var = predict2([0], x_train, y_train, kernel, output_noise = True)

print(y_mean)
[4.11090557]
print(y_var)
[0.00015]

今回は、カーネル行列\bm Kはきちんと正則行列となったようです。

実装する

以上の考察をもとに、ガウス過程回帰モデルをGaussianProcessRegressorクラスとして実装します。

class GaussianProcessRegressor:
    def __init__(self, kernel, sigma = 1e-2, output_noise = True):
        self.x_train = np.array([])
        self.y_train = np.array([])
        self.kernel = kernel
        self.sigma = sigma
        self.output_noise = output_noise

    def append(self, x_train, y_train):
        self.x_train = np.append(self.x_train, x_train)
        self.y_train = np.append(self.y_train, y_train)
    
    def kernel_matrix(self, xi, xj):
        xis, xjs = np.meshgrid(xj, xi)
        return self.kernel(xis, xjs)

    def predict(self, xs):
        K = self.kernel_matrix(self.x_train, self.x_train) \
            + self.sigma**2 * np.identity(len(self.x_train))
        K_inv = np.linalg.inv(K)

        y_mean = np.array([])
        y_var = np.array([])

        for x in xs:
            k_ = self.kernel_matrix(self.x_train, x)
            k__ = self.kernel_matrix(x, x)

            k_K = k_.T.dot(K_inv)
            y_mean = np.append(y_mean, k_K.dot(self.y_train))
            y_var = np.append(y_var, k__ - k_K.dot(k_) + self.sigma**2 * self.output_noise)
        return y_mean, y_var

訓練前

訓練する前の状態では、次のようになります。

GP = GaussianProcessRegressor(
    kernel = rbf([1.]),
    sigma = 0.01
    output_noise = True
)

x_min = 0
x_max = 10

x = np.linspace(x_min, x_max, 100)

y_mean, y_var = GP.predict(x)
y_std = np.sqrt(np.maximum(y_var, 0))

plt.figure(figsize=(12, 6))
plt.plot(x, y_mean, label='mean', color='#0066cc')
plt.fill_between(x, y_mean-y_std, y_mean+y_std, label='2σ credible region', color='#0066cc', alpha=0.3)
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.legend()
plt.show()

output

訓練例1

これまで利用してきたものと同じ関数fから訓練点X_\textrm{train}を取得し、それをもとに回帰してみます。

x = np.linspace(x_min, x_max, 100)
y = func(x)

x_train = np.linspace(x_min, x_max, num_init)
y_train = func(x_train)

GP = GaussianProcessRegressor(
    kernel = rbf([1.]),
    sigma = 0.01,
    output_noise = True
)

GP.append(x_train, y_train)

y_mean, y_var = GP.predict(x)
y_std = np.sqrt(np.maximum(y_var, 0))

plt.figure(figsize=(12, 6))
plt.plot(x, y, label='true', linestyle='dashed', color='black')
plt.plot(x, y_mean, label='mean', color='#0066cc')
plt.fill_between(x, y_mean-y_std, y_mean+y_std, label='2σ credible region', color='#0066cc', alpha=0.3)
plt.scatter(GP.x_train, GP.y_train, label='data', color='black', marker='x', s=200)
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.legend()
plt.show()

output

訓練例2

先程得た図では、等間隔で訓練点を得ているためにあまりおもしろい結果にはなりませんでした。今度は、ある程度偏りをもたせて訓練点を定め、回帰を行ってみましょう。

x = np.linspace(x_min, x_max, 100)
y = func(x)

x_train = np.array([0.5, 1.0, 1.5, 2.0, 2.5, 7.5, 8.0, 8.5, 9.0, 9.5])
y_train = func(x_train)

GP = GaussianProcessRegressor(
    kernel = rbf([1.]),
    sigma = 0.01,
    output_noise = True
)

GP.append(x_train, y_train)

y_mean, y_var = GP.predict(x)
y_std = np.sqrt(np.maximum(y_var, 0))

plt.figure(figsize=(12, 6))
plt.plot(x, y, label='true', linestyle='dashed', color='black')
plt.plot(x, y_mean, label='mean', color='#0066cc')
plt.fill_between(x, y_mean-y_std, y_mean+y_std, label='2σ credible region', color='#0066cc', alpha=0.3)
plt.scatter(GP.x_train, GP.y_train, label='data', color='black', marker='x', s=200)
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.legend()
plt.show()

output

訓練点が密集している範囲ではほぼ元の関数通りに回帰されていますが、真ん中付近は訓練点がほとんどないために信用区間が広がっています。

ここで、さらに真ん中にひとつ点を打ってみます。

GP.append(5, func(5))

y_mean, y_var = GP.predict(x)
y_std = np.sqrt(np.maximum(y_var, 0))

plt.figure(figsize=(12, 6))
plt.plot(x, y, label='true', linestyle='dashed', color='black')
plt.plot(x, y_mean, label='mean', color='#0066cc')
plt.fill_between(x, y_mean-y_std, y_mean+y_std, label='2σ credible region', color='#0066cc', alpha=0.3)
plt.scatter(GP.x_train, GP.y_train, label='data', color='black', marker='x', s=200)
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.legend()
plt.show()

output

相変わらず関数の形からズレてはいますが、新たに追加された訓練点付近については、先程よりも信用区間が狭くなっていることがわかると思います。このように、ガウス過程回帰モデルでは、もとの関数についての情報が少ない範囲では、信用区間が広くとられる傾向があり、これによって、どれくらい情報が足りていないのかを数値的に明示することができます。

追加で行うべきこと

本記事で実装したGaussianProcessRegressorクラスは簡易的なものです。まともな回帰分析のために用いるには不十分であり、追加で次のようなものの実装が必要です。

  • ハイパーパラメータチューニング - 非常に重要。単純な勾配降下法を利用したものだけでも、あったほうがよい。
  • 多変数関数への拡張 - 2変数関数f(x_1, x_2)への拡張から始めて、多変数関数f(\bm x)へ一般化する。
  • 計算コストの節約 - 精度を大きく落とすことなく計算コストを小さくしたい。補助変数法が有名。

参考文献

  1. Carl Edward Rasmussen and Christopher K. I. Williams: "Gaussian Process for Machine Learning," 2. Regression, https://gaussianprocess.org/gpml/chapters/, The MIT Press, 2006.
  2. 大羽成征, 持橋大地: 「ガウス過程と機械学習」, 機械学習プロフェッショナルシリーズ, 講談社, 2019.

Discussion

2022/8/21 追加

カーネル行列を計算する関数kernel_matrix()について変更。もとは

def kernel_matrix(xis, xjs, kernel):
    N1 = len(xis)
    N2 = len(xjs)
    return np.array(
        [kernel(xi, xj) for xi in xis for xj in xjs]
    ).reshape(N1, N2)

でしたが、これだと引数xisxjsともにリスト型でないと受け付けられなかったので、代わりに

def kernel_matrix(xi, xj, kernel):
    xis, xjs = np.meshgrid(xj, xi)
    return kernel(xis, xjs)

に変更しました。それに合わせて、predict()関数も若干変更しました。

ログインするとコメントできます