先回のカーネル法による正則化最小二乗法（１）で導出した
により、カーネル法を使った線形回帰を解いてみたいと思います。
式にあてはめて解くだけなので早速結果です。回帰のモデルは多項式としました。
ちゃんと近似できていますね。訓練データは正則化最小二乗法と同じにしたので、求められた曲線も同じとなりました。
結果は同じものになりましたが、カーネル法はカーネル法による正則化最小二乗法の最後で書いた通り、計算量に優位性があります。試しに多項式の次元数をいくつかふって、通常の正則化最小二乗法とカーネル法を用いた場合で計算時間がどれくらい変わるか試してみました。
日本語は文字化けしたので凡例は超適当な英語なんですけど、まあ伝わりますよね。計算時間はあくまでも僕の環境においてですが、ずいぶん違います。4万次元では100倍以上差が出ています。さらに、通常の正則化最小二乗法で5万次元まで増やすと、僕の環境ではメモリエラーになってしまいましたが、カーネル法を使えば計算可能でした。計算時間の違いと、そもそも計算できるかどうかっていう点でカーネル法の優位性が確認できました。
今回は多項式モデルを考え、数万次元の特徴ベクトルを持つ回帰も計算可能であるとわかりましたが、カーネル関数によっては数万次元どころか無限次元を扱うこともできます。参考：ガウスカーネル
今回のコードです。
# カーネル法による線形回帰 import matplotlib.pyplot as plt import numpy as np # カーネル関数は特徴ベクトルの内積 # 今回の特徴ベクトルはM次の多項式 def k(x0, x1): ret = 0 for i in range(M+1): ret += x0**i * x1**i return ret # 行列K（グラム行列） def gram_matrix(x): gram = np.empty([N, N]) for i in range(N): for j in range(N): gram[i, j] = k(x[i], x[j]) return gram # ランダムシードを固定 np.random.seed(0) # 多項式の最大べき乗数(x^0+...+x^M) M = 9 # 訓練データ数 N = 10 # 正則化係数λ（参考書に倣った値） lam = np.exp(-18) # 訓練データ x = np.linspace(0, 1, N) t = np.sin(2*np.pi*x.T) + np.random.normal(0, 0.2, N) # 係数aを求める a = np.linalg.inv(gram_matrix(x) + lam*np.eye(N)) @ t # 新たな入力x2に対する予測値yを求める x2 = np.linspace(0, 1, 100) y = np.empty(100) for i in range(x2.size): y[i] = np.vectorize(k)(x.T, np.repeat(x2[i], N)) @ a # 結果の表示 plt.xlim(0.0, 1.0) plt.ylim(-1.5, 1.5) plt.scatter(x, t) plt.plot(x2, y) plt.show()
