サポートベクトルマシンにおけるカーネルトリックについて理解しようとしたら理屈が多くて大変だった（その２：カーネル回帰、リッジ回帰等）

はじめに

　機械学習で用いる推論モデルを一通り理解しようとしています。今回は、カーネルトリックという手法についてまとめました。
　前回はざっくりとした概要と実装例をまとめています。
　

サポートベクトルマシンについて丁寧に理解しようとした（その１：MakeMoonsを例にして多項式/RBFカーネルを試した）
https://qiita.com/Fumio-eisan/items/f2553266a509d6f2d3a3

　今回のカーネルトリックという手法は私にとって非常に理解が大変な内容が多く苦労しました。ロジスティク回帰やニューラルネットワークを理解する以上に、数学的な知識が必要なように感じました。。
　なるべく順を追って理解が難しい点も踏まえて分かりやすくまとめていければと思います。

　概要は下記です。

• 非線形な関数を作りたい、ではどうしよう？
• カーネル関数を理解する
• ガウスカーネルを理解する

非線形な関数を作りたい、ではどうしよう？

　サポートベクトルマシンでは、非線形な境界について分類する（＝超平面分離）ことができる分類器です。
　下記左図のように分けたい境界がある場合、その境界に沿った関数を立式して区別します。
　さて、非線形な関数を作るテクニックは、（私が認識している範囲で）3つあります。一つはテイラー展開と呼ばれる手法で、$sin,cos$等三角関数を高次の関数で近似する方法です。次に、フーリエ級数と呼ばれる手法です。波動に関わるテーマで良く用いられる手法で、周期性のある波に対して、同じく周期性のある$sin,cos$の和で表すものです。
　そして最後がカーネル法と呼ばれるものです。これは、指数関数であるexpや多項式等の関数で表現する手法です。

　そして、このカーネル法と呼ばれる手法で非線形な回帰を行います。このとき、下記図のようにテストデータに合った関数を作ることを考えます。このときに見つけてほしい関数と似ても似つかない曲線が描かれる場合があります。これを、過学習と呼びます。対策として、正則化項を入れて過学習を防ぐことをします。

このような流れで最適な関数を作ります。
次から、もう少しかみ砕きつつ、実装をしながら理解していきたいと思います。

カーネル回帰を理解する

　それでは、カーネル回帰について進めます。先ほども軽く触れましたが、非線形な線を描くための手段と理解しています。似ている（と思っている）表現に、フーリエ級数があります。
　波動をsin,cos波の足し合わせで表現する考え方です。私自身は仕事において、ガス配管に流れる脈動を調べるためにこの手法を用いて解析をしていました。従い、それとカーネル回帰は似たロジックだなと思いました。
　

さて、カーネル回帰は$N$個のカーネル関数の和として表されます。この$N$個とは、テストデータのサンプル数になります。

f({\bf x}) = \sum_{i=1}^{N} \alpha_i k({\bf x}^{(i)}, {\bf x})

この$k({\bf x}^{(i)}, {\bf x})$で表されるところがカーネル関数と呼びます。このカーネル関数は何種類かあります。代表的なものにガウスカーネルと多項式カーネルがあります。以下のような式で表されます。

ガウスカーネル:
k({\bf x}, {\bf x}') = exp(-\beta \|{\bf x} - {\bf x}'\|^2)\\

多項式カーネル:
k({\bf x}, {\bf x}') =  ({\bf x} {\bf x}'+c)^p

ガウスカーネルをいじくりまわす

今回は、ガウスカーネルを使って表現したいと思います。ガウスカーネルは$β$の値が絶対値の高さ、$x'$がサンプル数の位置を示すと捉えて頂いて良いと思います。下記に$β$と'x''を変化させた際の関数を示します。

kernel.ipynb

import numpy as np

def kernel(xi, xj, beta=1):
    return np.exp(- beta * np.sum((xi - xj)**2))

X = np.arange(-1.5, 1.5, 0.05)
Y = np.zeros(len(X))
for i in range(len(X)):
    Y[i] = kernel(X[i], 0, 10)

X1 = np.arange(-1.5, 1.5, 0.05)
Y1 = np.zeros(len(X1))
for i in range(len(X1)):
    Y1[i] = kernel(X[i], 0, 1)

X2 = np.arange(-1.5, 1.5, 0.05)
Y2 = np.zeros(len(X2))
for i in range(len(X2)):
    Y2[i] = kernel(X[i], -1, 1)


plt.figure(figsize=(7, 4))
plt.plot(X2, Y2,label='beta=1,x=-1')
plt.plot(X1, Y1,label='beta=1')
plt.plot(X, Y,label='beta=10')
plt.legend()
plt.show()

$β$を大きくするとよりシャープな関数が描けることが分かります。また、$x'$の値を操作することで頂点となる位置がずれることが分かります。

　さて、次にこのカーネル関数の合成についてです。実際の非線形な関数を表現するうえでは、このガウスカーネルを足し合わせることが必要になります。
　下記プログラムでは、それぞれ生成されるガウスカーネルを足し合わせてプロットすることを行っています。

kernel.ipynb

def kernel(xi, xj, beta=1):
    return np.exp(- beta * np.sum((xi - xj)**2))

X = np.arange(-1.5, 1.5, 0.01)
Y = np.zeros(len(X))

centers = [0,-1,1]
weights = [1,-1,2]
for i in range(len(X)):
    for weight, center in zip(weights, centers):
        Y[i] += weight * kernel(X[i], center, 10)

plt.figure(figsize=(7, 4))
plt.plot(X, Y)
plt.show()

適当な数値を入れて計算させましたが、適当にぐちゃぐちゃした非線形な形を表現することができました。

テストデータに沿った最適化したカーネル関数を作る

　これで非線形な関数を得られることが分かりましたので、次に任意の非線形関数を作る手順について理解したいと思います。
　線形回帰の場合は、誤差関数の最小化により係数の最適化を行います。カーネル回帰においても考え方は一緒です。
　サンプル数$N$個のテストデータ（例：$(x^{(1)},y^{(1)})$）において、$alpha_{0～N}$の$N$個について最適化を行うことで非線形を表現可能な関数を作ることができます。
　

\hat{y} =
\left(\begin{matrix}
k({\bf x}^{(1)}, {\bf x}) & k({\bf x}^{(2)}, {\bf x}) & \cdots & k({\bf x}^{(N)}, {\bf x})
\end{matrix}\right)
\left(\begin{matrix}
\alpha_0 \\
\alpha_1 \\
\vdots \\
\alpha_N
\end{matrix}\right)

そして、

{\bf y} = 
\left(\begin{matrix}
y^{(1)} \\
y^{(2)} \\
\vdots \\
y^{(N)}
\end{matrix}\right)
,\;
K=\left(\begin{matrix}
k({\bf x}^{(1)}, {\bf x}^{(1)}) & k({\bf x}^{(2)}, {\bf x}^{(1)}) & \cdots & k({\bf x}^{(N)}, {\bf x}^{(1)}) \\
k({\bf x}^{(1)}, {\bf x}^{(2)}) & k({\bf x}^{(2)}, {\bf x}^{(2)}) & & \vdots \\
\vdots & & \ddots & \vdots \\
k({\bf x}^{(1)}, {\bf x}^{(N)}) & \cdots & \cdots & k({\bf x}^{(N)}, {\bf x}^{(N)})
\end{matrix}\right) 
,\;
{\bf \alpha}
=
\left(\begin{matrix}
\alpha^{(1)} \\
\alpha^{(2)} \\
\vdots \\
\alpha^{(N)}
\end{matrix}\right)

　誤差関数は$R(α)$で表すこととし、テストデータの値$y$と$f(x)$の差分の二乗としましょう。

\begin{align}
R({\bf \alpha}) &= \sum_{i=1}^{N}\left\{y^{(i)}- f({\bf x}^{(i)})\right\}^2
=\sum_{i=1}^{N}\left\{y^{(i)}- \sum_{i=1}^{N} \alpha_i k({\bf x}^{(i)}, {\bf x})\right\}^2 \\
&=(\begin{array} yy_1 - \alpha_1k(\bf x^{(1)}-\bf x)&y_2 - \alpha_2k(\bf x^{(2)}-\bf x)  &\cdots&y_n - \alpha_Nk(\bf x^{(N)}-\bf x) \end{array})
\begin{pmatrix}
y_1 - \alpha_1k(\bf x^{(1)}-\bf x)\\
y_2 - \alpha_2k(\bf x^{(2)}-\bf x) \\
\vdots \\
y_n - \alpha_Nk(\bf x^{(N)}-\bf x)\\ 
\end{pmatrix}\\

&=({\bf y}- K {\bf \alpha})^{\mathrm{T}}({\bf y}- K {\bf \alpha})
\end{align}

にように誤差関数を表すことができます。また、$\alpha$に関しては、

\bf y = \bf K\bf α \\
∴\bf α = \bf K^{-1} \bf y

として、カーネル関数の逆行列を用いることで求めることができます。

さて、適当なデータセットとして$cos$関数で描かれる曲線を再現させたいと思います。

kernel.ipynb

import pandas as pd

def wave_dataset(size, xlim=[0, 2], scale=None):
    x = np.random.uniform(xlim[0], xlim[1], size)
    y = np.cos(x)
    if scale is not None:
        noize = np.random.normal(0, scale, size)
        y = y + noize
    df = pd.DataFrame()
    df['x'] = x
    df['y'] = y
    return df

data = wave_dataset(20, xlim=[-3.14, 3.14], scale=0.05)

plt.figure(figsize=(7, 4))
plt.scatter(data['x'].values, data['y'].values)
plt.ylim((-1.5, 1.5))
plt.show()

np.random.uniform()メソッドによってランダムな値を生成することができます。この値を$cos$関数により得られた値に足して多少バラバラな値を生成しました。

次に、カーネル関数に値を入れていきます。そして、$\alpha$を逆行列から計算します。

kernel.ipynb

from itertools import combinations_with_replacement

X = data['x'].values
Y = data['y'].values

# ハイパーパラメータ
beta = 1

# データ数
N = X.shape[0]

# グラム行列の計算
K = np.zeros((N, N))
for i, j in combinations_with_replacement(range(N), 2):
    K[i][j] = kernel(X[i], X[j])
    K[j][i] = K[i][j]

# 重みを計算
alpha = np.linalg.inv(K).dot(Y)

これで準備が整いました。最後に、データセットを用いてカーネル回帰をして関数を生成・表記してみます。

kernel.ipynb

# カーネル回帰
def kernel_predict(X, x, alpha, beta):
    Y = 0
    for i in range(len(X)):
        Y += alpha[i] * kernel(X[i], x, beta)
    return Y

# 回帰によって結果を予測        
X_axis = np.arange(-3.14, 3.14, 0.01)
Y_predict = np.zeros(len(X_axis))
for i in range(len(X_axis)):
    Y_predict[i] = kernel_predict(X, X_axis[i], alpha, beta)

# 結果を描画
plt.figure(figsize=(7, 4))

## 観測データ
plt.scatter(data['x'].values, data['y'].values)
## 真の関数
plt.plot(X_axis, np.cos(X_axis), c='red')
## 予測した関数
plt.plot(X_axis, Y_predict)

plt.ylim((-1.5, 1.5))
plt.show()

青がカーネル回帰により生成した関数です。確かにテストデータの点に近いところを通過していますが、$cos$関数とは異なる形になってしまいました。
　これは、モデルの過学習が問題です。実は$α$の値は下記のようになっています。

$10^{9～13}$のオーダーの非常に大きな値が入っていることが分かります。これでは、得られる最大値もそれだけ大きな値を取ってしまいます。

対策として、正則化項を入れる手続きを行うことで解決します。

リッジ回帰を理解する

　さて、正則化項とは、先ほどの二乗誤差に対して、正則化項と呼ばれる項を付け加えて最適化することをいいます。

\begin{align}
R({\bf \alpha})
&=({\bf y}- K {\bf \alpha})^{\mathrm{T}}({\bf y}- K {\bf \alpha}) + \lambda {\bf \alpha}^{\mathrm{T}} K {\bf \alpha}
\end{align}

$\lambda {\bf \alpha}^{\mathrm{T}} K {\bf \alpha}$を加えて損失関数の最小化を行うことをリッジ回帰と呼びます。この項のことをリッジ回帰の正則化項と呼びます。
　基本的にこの損失関数の値を最小化させたいため、必然的に$α$も小さくする必要があります。
　この条件のもとで$\alpha$を求めると以下のように表すことができます。

{\bf \alpha} = (K+\lambda I_N)^{-1}{\bf y}

この考え方を用いて先ほどの$\alpha$を求めなおして、再度プロットさせてみたいと思います。

kernel.ipynb

# 正則化項の係数
lam = 0.5
alpha_r = np.linalg.inv(K + lam * np.eye(K.shape[0])).dot(Y)

# 回帰によって結果を予測
Y_predict_r = np.zeros(len(X_axis))
for i in range(len(X_axis)):
    Y_predict_r[i] = kernel_predict(X, X_axis[i], alpha_r, beta)

# 結果を描画
plt.figure(figsize=(6.472, 4))

## 観測データ
plt.scatter(data['x'].values, data['y'].values)
## 真の関数
plt.plot(X_axis, np.cos(X_axis), c='red')
## 予測した関数
plt.plot(X_axis, Y_predict_r)

plt.ylim((-1.5, 1.5))
plt.show()

できました。凡そ再現できたといっていいかと思います（$x=3$の辺りはテストデータが少ないこともあり変曲点となっていますが。。）。

おわりに

　駆け足になったかもしれませんが、サポートベクトルマシンにおける非線形関数を作る手続きをまとめました。数学的に理解しなければならないことが多く、大変でした。
　しかし、一方でその考え方が分かることでこのアルゴリズムの魅力を知れたような気がします。

多大に参考にさせて頂いた記事がこちらです。

線形な手法とカーネル法（回帰分析）
https://qiita.com/wsuzume/items/09a59036c8944fd563ff

　
プログラム全文はこちらです。
https://github.com/Fumio-eisan/SVM_20200417