カーネルリッジ回帰（Kernel Ridge Regression, KRR）は、リッジ回帰をカーネル法によって非線形に拡張した手法です。その理論的基盤は再生核ヒルベルト空間（RKHS）にあり、リプレゼンター定理によって「最適解はデータ点のカーネル関数の線形結合で表される」ことが保証されます。
さらに、カーネルリッジ回帰はガウス過程回帰と深い関係があります。ガウス過程の事前分布が RKHS のノルムによる正則化に対応するという意味で、両者は同一の予測関数を与えます。本記事では、リッジ回帰の復習から出発し、リプレゼンター定理の証明、カーネルリッジ回帰の閉形式解の導出、ガウス過程回帰との等価性の議論、計算量の比較までを丁寧に解説し、Python での実装と比較を行います。
本記事の内容
- リッジ回帰の復習と特徴空間への拡張
- リプレゼンター定理の主張と証明
- カーネルリッジ回帰の閉形式解 $\bm{\alpha}^* = (\bm{K} + \lambda\bm{I})^{-1}\bm{y}$
- ガウス過程回帰との等価性
- 計算量の比較（$O(N^3)$ vs $O(D^3)$）
- Python での各種カーネルによる実装とガウス過程回帰との比較
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
リッジ回帰の復習
リッジ回帰の定式化
訓練データ $\{(\bm{x}_i, y_i)\}_{i=1}^{N}$（$\bm{x}_i \in \mathbb{R}^d$, $y_i \in \mathbb{R}$）に対して、リッジ回帰は以下の正則化付き最小二乗問題を解きます。
$$ \min_{\bm{w}} \frac{1}{2}\sum_{i=1}^{N}(y_i – \bm{w}^T\bm{x}_i)^2 + \frac{\lambda}{2}\|\bm{w}\|^2 $$
行列形式では、$\bm{X} = [\bm{x}_1, \dots, \bm{x}_N]^T \in \mathbb{R}^{N \times d}$, $\bm{y} = (y_1, \dots, y_N)^T$ として
$$ \min_{\bm{w}} \frac{1}{2}\|\bm{y} – \bm{X}\bm{w}\|^2 + \frac{\lambda}{2}\|\bm{w}\|^2 $$
閉形式解の導出
目的関数を $J(\bm{w})$ とおき、$\bm{w}$ で微分します。
$$ \nabla_{\bm{w}} J = -\bm{X}^T(\bm{y} – \bm{X}\bm{w}) + \lambda\bm{w} = \bm{0} $$
整理すると
$$ (\bm{X}^T\bm{X} + \lambda\bm{I}_d)\bm{w} = \bm{X}^T\bm{y} $$
$$ \bm{w}^* = (\bm{X}^T\bm{X} + \lambda\bm{I}_d)^{-1}\bm{X}^T\bm{y} $$
$\lambda > 0$ により $\bm{X}^T\bm{X} + \lambda\bm{I}_d$ は常に正則であり、逆行列が確実に存在します。計算量は $O(d^3 + Nd^2)$ で、$d \times d$ 行列の逆行列計算が支配的です。
予測
新しいデータ点 $\bm{x}_*$ での予測は
$$ f(\bm{x}_*) = \bm{w}^{*T}\bm{x}_* = \bm{y}^T\bm{X}(\bm{X}^T\bm{X} + \lambda\bm{I}_d)^{-1}\bm{x}_* $$
Woodbury の恒等式による双対形式
ここで重要な行列恒等式を用います。Woodbury の恒等式（push-through identity の一種）により
$$ \bm{X}(\bm{X}^T\bm{X} + \lambda\bm{I}_d)^{-1}\bm{X}^T = (\bm{X}\bm{X}^T + \lambda\bm{I}_N)^{-1}\bm{X}\bm{X}^T $$
これを用いると、予測関数は
$$ f(\bm{x}_*) = \bm{y}^T(\bm{X}\bm{X}^T + \lambda\bm{I}_N)^{-1}\bm{X}\bm{X}^T \cdot \bm{X}(\bm{X}^T\bm{X}+\lambda\bm{I}_d)^{-1}\bm{x}_* $$
ではなく、より直接的に以下のように書けます。$\bm{w}^* = \bm{X}^T\bm{\alpha}$ と分解することで
$$ f(\bm{x}_*) = \bm{k}_*^T(\bm{K} + \lambda\bm{I}_N)^{-1}\bm{y} $$
ここで $\bm{K} = \bm{X}\bm{X}^T$（$K_{ij} = \bm{x}_i^T\bm{x}_j$）はグラム行列、$\bm{k}_* = \bm{X}\bm{x}_*$（$(\bm{k}_*)_i = \bm{x}_i^T\bm{x}_*$）です。
この形式の計算量は $O(N^3)$（$N \times N$ 行列の逆行列）です。$d \gg N$ のとき、あるいは $d = \infty$（カーネル法）のとき、この双対形式が有利になります。
特徴空間でのリッジ回帰
特徴写像の導入
特徴写像 $\phi: \mathbb{R}^d \to \mathcal{H}$ を導入し、特徴空間でのリッジ回帰を考えます。
$$ \min_{\bm{w}} \frac{1}{2}\sum_{i=1}^{N}(y_i – \bm{w}^T\phi(\bm{x}_i))^2 + \frac{\lambda}{2}\|\bm{w}\|^2 $$
特徴行列 $\bm{\Phi} = [\phi(\bm{x}_1), \dots, \phi(\bm{x}_N)]^T \in \mathbb{R}^{N \times D}$（$D$ は特徴空間の次元）を用いると
$$ \min_{\bm{w}} \frac{1}{2}\|\bm{y} – \bm{\Phi}\bm{w}\|^2 + \frac{\lambda}{2}\|\bm{w}\|^2 $$
主問題の解は
$$ \bm{w}^* = (\bm{\Phi}^T\bm{\Phi} + \lambda\bm{I}_D)^{-1}\bm{\Phi}^T\bm{y} $$
しかし $D$ が非常に大きい（あるいは無限の）場合、$D \times D$ 行列の逆行列計算は不可能です。
ここでカーネルトリックの出番です。双対形式を用いれば、$N \times N$ のカーネル行列 $K_{ij} = k(\bm{x}_i, \bm{x}_j) = \phi(\bm{x}_i)^T\phi(\bm{x}_j)$ だけで計算できます。しかし、なぜ最適解が双対形式で書けるのでしょうか？ その理論的根拠がリプレゼンター定理です。
リプレゼンター定理
定理の主張
リプレゼンター定理（Representer Theorem）: $k$ を正定値カーネル、$\mathcal{H}_k$ を対応する RKHS とする。正則化付き経験リスク最小化問題
$$ \min_{f \in \mathcal{H}_k} \sum_{i=1}^{N} L(y_i, f(\bm{x}_i)) + \lambda \|f\|_{\mathcal{H}_k}^2 $$
の最適解 $f^*$ は、以下の形で表される。
$$ f^*(\bm{x}) = \sum_{i=1}^{N} \alpha_i^* k(\bm{x}_i, \bm{x}) $$
すなわち、最適解はデータ点におけるカーネル関数の線形結合で表されます。$L$ は任意の損失関数（微分可能性すら不要）であり、$\lambda > 0$ は正則化パラメータです。
証明
RKHS $\mathcal{H}_k$ の任意の関数 $f$ は、データ点に関連する成分とそれに直交する成分に分解できます。
$\mathcal{V} = \text{span}\{k(\bm{x}_1, \cdot), \dots, k(\bm{x}_N, \cdot)\}$ を、データ点でのカーネル関数が張る $\mathcal{H}_k$ の部分空間とします。
任意の $f \in \mathcal{H}_k$ は、直交分解
$$ f = f_{\parallel} + f_{\perp} $$
で表せます。ここで $f_{\parallel} \in \mathcal{V}$、$f_{\perp} \in \mathcal{V}^{\perp}$ です。
Step 1: $f_{\parallel}$ は以下の形で書けます。
$$ f_{\parallel}(\bm{x}) = \sum_{i=1}^{N} \alpha_i k(\bm{x}_i, \bm{x}) $$
Step 2: 再生性 $f(\bm{x}_j) = \langle f, k(\bm{x}_j, \cdot) \rangle_{\mathcal{H}_k}$ を用いると、データ点での関数値は
$$ f(\bm{x}_j) = \langle f_{\parallel} + f_{\perp}, k(\bm{x}_j, \cdot) \rangle_{\mathcal{H}_k} = \langle f_{\parallel}, k(\bm{x}_j, \cdot) \rangle_{\mathcal{H}_k} + \langle f_{\perp}, k(\bm{x}_j, \cdot) \rangle_{\mathcal{H}_k} $$
$k(\bm{x}_j, \cdot) \in \mathcal{V}$ であり、$f_{\perp} \in \mathcal{V}^{\perp}$ なので
$$ \langle f_{\perp}, k(\bm{x}_j, \cdot) \rangle_{\mathcal{H}_k} = 0 $$
よって
$$ f(\bm{x}_j) = f_{\parallel}(\bm{x}_j), \quad j = 1, \dots, N $$
つまり、データ点での関数値は $f_{\perp}$ に依存しません。
Step 3: 損失関数 $L(y_i, f(\bm{x}_i))$ はデータ点での関数値のみに依存するので
$$ \sum_{i=1}^{N} L(y_i, f(\bm{x}_i)) = \sum_{i=1}^{N} L(y_i, f_{\parallel}(\bm{x}_i)) $$
Step 4: 正則化項については、直交分解の性質から
$$ \|f\|_{\mathcal{H}_k}^2 = \|f_{\parallel}\|_{\mathcal{H}_k}^2 + \|f_{\perp}\|_{\mathcal{H}_k}^2 \geq \|f_{\parallel}\|_{\mathcal{H}_k}^2 $$
等号は $f_{\perp} = 0$ のとき、つまり $f = f_{\parallel}$ のときにのみ成立します。
Step 5: 以上より、目的関数は
$$ \sum_{i=1}^{N} L(y_i, f(\bm{x}_i)) + \lambda\|f\|_{\mathcal{H}_k}^2 = \sum_{i=1}^{N} L(y_i, f_{\parallel}(\bm{x}_i)) + \lambda(\|f_{\parallel}\|_{\mathcal{H}_k}^2 + \|f_{\perp}\|_{\mathcal{H}_k}^2) $$
$\|f_{\perp}\|_{\mathcal{H}_k}^2 \geq 0$ であり、第1項は $f_{\perp}$ に依存しないため、$f_{\perp} = 0$ とすることで目的関数は減少（少なくとも非増加）します。
したがって、最適解は $f^* = f_{\parallel}^* \in \mathcal{V}$ の形、すなわち
$$ f^*(\bm{x}) = \sum_{i=1}^{N} \alpha_i^* k(\bm{x}_i, \bm{x}) $$
で表されます。$\blacksquare$
リプレゼンター定理の意義
リプレゼンター定理により、無限次元空間 $\mathcal{H}_k$ での最適化問題が、有限（$N$ 個）のパラメータ $\bm{\alpha} = (\alpha_1, \dots, \alpha_N)^T$ の最適化問題に帰着されます。これは計算上極めて重要です。
カーネルリッジ回帰の閉形式解
導出
リプレゼンター定理より、最適解は $f(\bm{x}) = \sum_{i=1}^{N} \alpha_i k(\bm{x}_i, \bm{x})$ の形です。これをリッジ回帰の目的関数に代入します。
データ点での予測値は
$$ f(\bm{x}_j) = \sum_{i=1}^{N} \alpha_i k(\bm{x}_i, \bm{x}_j) = (\bm{K}\bm{\alpha})_j $$
ベクトル形式で $\bm{f} = \bm{K}\bm{\alpha}$ です。
RKHS ノルムは、$f = \sum_i \alpha_i k(\bm{x}_i, \cdot)$ のとき
$$ \|f\|_{\mathcal{H}_k}^2 = \left\langle \sum_i \alpha_i k(\bm{x}_i, \cdot), \sum_j \alpha_j k(\bm{x}_j, \cdot) \right\rangle_{\mathcal{H}_k} = \sum_i \sum_j \alpha_i \alpha_j k(\bm{x}_i, \bm{x}_j) = \bm{\alpha}^T \bm{K} \bm{\alpha} $$
ここで再生性 $\langle k(\bm{x}_i, \cdot), k(\bm{x}_j, \cdot) \rangle_{\mathcal{H}_k} = k(\bm{x}_i, \bm{x}_j)$ を用いました。
目的関数は
$$ J(\bm{\alpha}) = \frac{1}{2}\|\bm{y} – \bm{K}\bm{\alpha}\|^2 + \frac{\lambda}{2}\bm{\alpha}^T\bm{K}\bm{\alpha} $$
$\bm{\alpha}$ で微分して 0 とおきます。
$$ \nabla_{\bm{\alpha}} J = -\bm{K}(\bm{y} – \bm{K}\bm{\alpha}) + \lambda\bm{K}\bm{\alpha} = \bm{0} $$
$$ -\bm{K}\bm{y} + \bm{K}^2\bm{\alpha} + \lambda\bm{K}\bm{\alpha} = \bm{0} $$
$$ \bm{K}(\bm{K} + \lambda\bm{I})\bm{\alpha} = \bm{K}\bm{y} $$
$\bm{K}$ が正定値の場合は両辺から $\bm{K}$ を消去できます。また、$\bm{K}$ が半正定値の場合でも $\bm{K} + \lambda\bm{I}$（$\lambda > 0$）は正定値なので
$$ \bm{\alpha}^* = (\bm{K} + \lambda\bm{I})^{-1}\bm{y} $$
が閉形式解となります。
新しいデータ点での予測
$$ f(\bm{x}_*) = \sum_{i=1}^{N} \alpha_i^* k(\bm{x}_i, \bm{x}_*) = \bm{k}_*^T \bm{\alpha}^* = \bm{k}_*^T (\bm{K} + \lambda\bm{I})^{-1}\bm{y} $$
ここで $\bm{k}_* = (k(\bm{x}_1, \bm{x}_*), \dots, k(\bm{x}_N, \bm{x}_*))^T$ です。
この式は特徴写像 $\phi$ を一切使わず、カーネル関数 $k$ のみで計算が完結しています。
ガウス過程回帰との等価性
ガウス過程回帰の復習
ガウス過程回帰（GPR）では、関数に事前分布としてガウス過程を仮定します。
$$ f \sim \mathcal{GP}(0, k(\bm{x}, \bm{x}’)) $$
観測モデル $y_i = f(\bm{x}_i) + \varepsilon_i$、$\varepsilon_i \sim \mathcal{N}(0, \sigma_n^2)$ のもとで、事後平均は
$$ \bar{f}(\bm{x}_*) = \bm{k}_*^T(\bm{K} + \sigma_n^2 \bm{I})^{-1}\bm{y} $$
等価性
カーネルリッジ回帰の予測とガウス過程回帰の事後平均を比較すると
| カーネルリッジ回帰 | ガウス過程回帰 | |
|---|---|---|
| 予測関数 | $\bm{k}_*^T(\bm{K} + \lambda\bm{I})^{-1}\bm{y}$ | $\bm{k}_*^T(\bm{K} + \sigma_n^2\bm{I})^{-1}\bm{y}$ |
| 正則化 / ノイズ | $\lambda$ | $\sigma_n^2$ |
$\lambda = \sigma_n^2$ とおけば、カーネルリッジ回帰の予測関数はガウス過程回帰の事後平均と完全に一致します。
理論的背景
この等価性はより深い理論的背景を持ちます。
RKHS の観点: カーネルリッジ回帰は、RKHS ノルム $\|f\|_{\mathcal{H}_k}$ による正則化を通じて、「滑らかな」関数を選好します。
ベイズの観点: ガウス過程では、事前分布のカーネル $k$ が関数の「滑らかさ」を定義し、事後分布は尤度（データへのフィット）と事前分布（滑らかさ）のバランスで決まります。
両者の関係は次のように整理できます。
$$ \underbrace{\|f\|_{\mathcal{H}_k}^2}_{\text{RKHS正則化}} \longleftrightarrow \underbrace{-2\log p(f)}_{\text{ガウス過程事前分布の対数}} $$
RKHS ノルムが小さい関数は、ガウス過程の事前分布のもとで高い確率を持つ関数に対応します。
違い：不確実性の定量化
ただし、両者には重要な違いもあります。
- カーネルリッジ回帰: 点予測のみ（$f(\bm{x}_*)$ の値）
- ガウス過程回帰: 予測分布を返す（平均 $\bar{f}(\bm{x}_*)$ に加えて分散 $\text{Var}(f(\bm{x}_*))$ も得られる）
ガウス過程回帰の事後分散は
$$ \text{Var}(f(\bm{x}_*)) = k(\bm{x}_*, \bm{x}_*) – \bm{k}_*^T(\bm{K} + \sigma_n^2\bm{I})^{-1}\bm{k}_* $$
であり、予測の不確実性を定量化できます。これはカーネルリッジ回帰の枠組みからは直接得られません。
計算量の比較
主問題 vs 双対問題
| 主問題（特徴空間） | 双対問題（カーネル法） | |
|---|---|---|
| 逆行列 | $(\bm{\Phi}^T\bm{\Phi} + \lambda\bm{I}_D)^{-1}$: $D \times D$ | $(\bm{K} + \lambda\bm{I}_N)^{-1}$: $N \times N$ |
| 訓練 | $O(ND^2 + D^3)$ | $O(N^2D_k + N^3)$ |
| 予測 | $O(D)$ per sample | $O(N D_k)$ per sample |
| 有利な場合 | $D \ll N$ | $N \ll D$ or $D = \infty$ |
ここで $D$ は特徴空間の次元、$D_k$ はカーネル計算に必要な元の入力次元、$N$ はデータ数です。
大規模データへの対応
$N$ が大きい場合（例えば $N > 10000$）、$N \times N$ のカーネル行列の逆行列計算 $O(N^3)$ は現実的ではありません。この場合、以下の近似手法が用いられます。
- ナイストレム近似: ランダムに選んだ $M \ll N$ 個のランドマーク点でカーネル行列を低ランク近似
- ランダム特徴: RFF（Random Fourier Features）により $\phi(\bm{x}) \in \mathbb{R}^{D_{RF}}$ を明示的に構成し、主問題として解く
- 疎近似: 誘導点（inducing points）を用いたガウス過程の疎近似
Python での実装と比較
カーネルリッジ回帰の実装
import numpy as np
import matplotlib.pyplot as plt
def rbf_kernel(X1, X2, length_scale=1.0):
"""RBFカーネル"""
sq_dists = np.sum(X1**2, axis=1, keepdims=True) - 2 * X1 @ X2.T + np.sum(X2**2, axis=1, keepdims=True).T
return np.exp(-0.5 * sq_dists / length_scale**2)
def matern32_kernel(X1, X2, length_scale=1.0):
"""マテルン3/2カーネル"""
dists = np.sqrt(np.sum(X1**2, axis=1, keepdims=True) - 2 * X1 @ X2.T +
np.sum(X2**2, axis=1, keepdims=True).T + 1e-12)
r = np.sqrt(3) * dists / length_scale
return (1 + r) * np.exp(-r)
def periodic_kernel(X1, X2, length_scale=1.0, period=2*np.pi):
"""周期カーネル"""
dists = np.sqrt(np.sum((X1[:, None, :] - X2[None, :, :])**2, axis=2) + 1e-12)
return np.exp(-2 * np.sin(np.pi * dists / period)**2 / length_scale**2)
class KernelRidgeRegression:
"""カーネルリッジ回帰"""
def __init__(self, kernel_func, lam=1.0, **kernel_params):
self.kernel_func = kernel_func
self.lam = lam
self.kernel_params = kernel_params
self.alpha = None
self.X_train = None
def fit(self, X, y):
"""学習: α = (K + λI)^{-1} y"""
self.X_train = X.copy()
K = self.kernel_func(X, X, **self.kernel_params)
self.alpha = np.linalg.solve(K + self.lam * np.eye(len(X)), y)
def predict(self, X):
"""予測: f(x*) = k*^T α"""
k_star = self.kernel_func(X, self.X_train, **self.kernel_params)
return k_star @ self.alpha
class GaussianProcessRegression:
"""ガウス過程回帰"""
def __init__(self, kernel_func, sigma_n=0.1, **kernel_params):
self.kernel_func = kernel_func
self.sigma_n = sigma_n
self.kernel_params = kernel_params
self.X_train = None
self.y_train = None
self.K_inv = None
def fit(self, X, y):
"""学習"""
self.X_train = X.copy()
self.y_train = y.copy()
K = self.kernel_func(X, X, **self.kernel_params)
self.K_inv = np.linalg.inv(K + self.sigma_n**2 * np.eye(len(X)))
def predict(self, X, return_std=False):
"""予測（事後平均と事後分散）"""
k_star = self.kernel_func(X, self.X_train, **self.kernel_params)
mean = k_star @ self.K_inv @ self.y_train
if return_std:
k_ss = self.kernel_func(X, X, **self.kernel_params)
var = np.diag(k_ss - k_star @ self.K_inv @ k_star.T)
var = np.maximum(var, 0) # 数値誤差対策
return mean, np.sqrt(var)
return mean
各種カーネルでのカーネルリッジ回帰
# 合成データの生成
np.random.seed(42)
N = 50
X_train = np.sort(np.random.uniform(-5, 5, N)).reshape(-1, 1)
y_true_func = lambda x: np.sin(x) + 0.3 * np.sin(3 * x)
y_train = y_true_func(X_train.ravel()) + 0.2 * np.random.randn(N)
X_test = np.linspace(-6, 6, 300).reshape(-1, 1)
y_test_true = y_true_func(X_test.ravel())
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
kernel_configs = [
(rbf_kernel, {'length_scale': 1.0}, 'RBF Kernel ($\\ell=1.0$)'),
(rbf_kernel, {'length_scale': 0.3}, 'RBF Kernel ($\\ell=0.3$)'),
(matern32_kernel, {'length_scale': 1.0}, 'Matern 3/2 Kernel ($\\ell=1.0$)'),
]
for ax, (kernel_func, params, title) in zip(axes, kernel_configs):
krr = KernelRidgeRegression(kernel_func, lam=0.01, **params)
krr.fit(X_train, y_train)
y_pred = krr.predict(X_test)
ax.scatter(X_train, y_train, c='gray', s=20, alpha=0.7, label='Training data')
ax.plot(X_test, y_test_true, 'g--', linewidth=2, label='True function')
ax.plot(X_test, y_pred, 'r-', linewidth=2, label='KRR prediction')
ax.set_title(title, fontsize=13)
ax.legend(fontsize=9)
ax.set_ylim(-2.5, 2.5)
ax.grid(True, alpha=0.3)
plt.suptitle('Kernel Ridge Regression with Different Kernels', fontsize=15, y=1.02)
plt.tight_layout()
plt.savefig('krr_different_kernels.png', dpi=150, bbox_inches='tight')
plt.show()
カーネルリッジ回帰 vs ガウス過程回帰の比較
# KRRとGPRの比較（λ = σ_n^2 で予測が一致することを確認）
sigma_n = 0.2
lam = sigma_n ** 2 # λ = σ_n^2
krr = KernelRidgeRegression(rbf_kernel, lam=lam, length_scale=1.0)
krr.fit(X_train, y_train)
y_pred_krr = krr.predict(X_test)
gpr = GaussianProcessRegression(rbf_kernel, sigma_n=sigma_n, length_scale=1.0)
gpr.fit(X_train, y_train)
y_pred_gpr, y_std_gpr = gpr.predict(X_test, return_std=True)
fig, axes = plt.subplots(1, 2, figsize=(16, 6))
# KRR
axes[0].scatter(X_train, y_train, c='gray', s=20, alpha=0.7, label='Training data')
axes[0].plot(X_test, y_test_true, 'g--', linewidth=2, label='True function')
axes[0].plot(X_test, y_pred_krr, 'r-', linewidth=2, label='KRR prediction')
axes[0].set_title(f'Kernel Ridge Regression ($\\lambda={lam}$)', fontsize=13)
axes[0].legend(fontsize=10)
axes[0].set_ylim(-2.5, 2.5)
axes[0].grid(True, alpha=0.3)
# GPR（不確実性付き）
axes[1].scatter(X_train, y_train, c='gray', s=20, alpha=0.7, label='Training data')
axes[1].plot(X_test, y_test_true, 'g--', linewidth=2, label='True function')
axes[1].plot(X_test, y_pred_gpr, 'b-', linewidth=2, label='GPR mean')
axes[1].fill_between(X_test.ravel(),
y_pred_gpr - 2 * y_std_gpr,
y_pred_gpr + 2 * y_std_gpr,
alpha=0.2, color='blue', label='$\\pm 2\\sigma$')
axes[1].set_title(f'Gaussian Process Regression ($\\sigma_n={sigma_n}$)', fontsize=13)
axes[1].legend(fontsize=10)
axes[1].set_ylim(-2.5, 2.5)
axes[1].grid(True, alpha=0.3)
plt.suptitle('KRR vs GPR: Same Prediction, Different Information', fontsize=15, y=1.02)
plt.tight_layout()
plt.savefig('krr_vs_gpr.png', dpi=150, bbox_inches='tight')
plt.show()
# 予測値の差を確認
diff = np.max(np.abs(y_pred_krr - y_pred_gpr))
print(f"KRRとGPRの予測の最大差: {diff:.2e}")
print(f"→ λ = σ_n^2 で予測は実質的に一致します")
正則化パラメータの影響
# 正則化パラメータ λ の影響
fig, axes = plt.subplots(2, 3, figsize=(16, 10))
lambdas = [1e-6, 1e-4, 0.01, 0.1, 1.0, 10.0]
for idx, lam in enumerate(lambdas):
ax = axes[idx // 3, idx % 3]
krr = KernelRidgeRegression(rbf_kernel, lam=lam, length_scale=1.0)
krr.fit(X_train, y_train)
y_pred = krr.predict(X_test)
ax.scatter(X_train, y_train, c='gray', s=15, alpha=0.7)
ax.plot(X_test, y_test_true, 'g--', linewidth=1.5, label='True')
ax.plot(X_test, y_pred, 'r-', linewidth=2, label='KRR')
ax.set_title(f'$\\lambda = {lam}$', fontsize=13)
ax.legend(fontsize=9)
ax.set_ylim(-3, 3)
ax.grid(True, alpha=0.3)
plt.suptitle('Effect of Regularization $\\lambda$ on KRR (RBF Kernel, $\\ell=1.0$)', fontsize=15, y=1.02)
plt.tight_layout()
plt.savefig('krr_lambda_effect.png', dpi=150, bbox_inches='tight')
plt.show()
リプレゼンター定理の確認
# リプレゼンター定理の確認: f*(x) = Σ α_i k(x_i, x) で本当に表現されているか
krr = KernelRidgeRegression(rbf_kernel, lam=0.04, length_scale=1.0)
krr.fit(X_train, y_train)
fig, axes = plt.subplots(1, 2, figsize=(16, 5))
# 各基底関数 α_i * k(x_i, ·) の可視化
ax = axes[0]
for i in range(0, N, 5): # 10個ごとに表示
basis = krr.alpha[i] * rbf_kernel(X_test, X_train[i:i+1], length_scale=1.0).ravel()
ax.plot(X_test, basis, alpha=0.4, linewidth=1)
ax.set_title('Individual basis functions $\\alpha_i k(x_i, \\cdot)$', fontsize=13)
ax.set_xlabel('$x$', fontsize=12)
ax.grid(True, alpha=0.3)
# 合計と予測の比較
ax = axes[1]
y_pred = krr.predict(X_test)
# 基底関数の合計を手動計算
y_manual = np.zeros(len(X_test))
for i in range(N):
y_manual += krr.alpha[i] * rbf_kernel(X_test, X_train[i:i+1], length_scale=1.0).ravel()
ax.plot(X_test, y_pred, 'r-', linewidth=2, label='KRR predict()')
ax.plot(X_test, y_manual, 'b--', linewidth=2, label='$\\sum_i \\alpha_i k(x_i, \\cdot)$')
ax.scatter(X_train, y_train, c='gray', s=20, alpha=0.5)
ax.set_title('Representer Theorem Verification', fontsize=13)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('representer_theorem.png', dpi=150, bbox_inches='tight')
plt.show()
# 数値的な一致を確認
print(f"predict()と手動計算の最大差: {np.max(np.abs(y_pred - y_manual)):.2e}")
計算量のスケーリング実験
# 計算量の比較（N vs 計算時間）
import time
Ns = [50, 100, 200, 500, 1000, 2000]
times_krr = []
times_gpr = []
for n in Ns:
X_temp = np.random.randn(n, 5) # 5次元データ
y_temp = np.random.randn(n)
# KRR
start = time.time()
krr = KernelRidgeRegression(rbf_kernel, lam=0.1, length_scale=1.0)
krr.fit(X_temp, y_temp)
times_krr.append(time.time() - start)
# GPR
start = time.time()
gpr = GaussianProcessRegression(rbf_kernel, sigma_n=np.sqrt(0.1), length_scale=1.0)
gpr.fit(X_temp, y_temp)
times_gpr.append(time.time() - start)
fig, ax = plt.subplots(figsize=(10, 6))
ax.loglog(Ns, times_krr, 'ro-', linewidth=2, markersize=8, label='KRR (solve)')
ax.loglog(Ns, times_gpr, 'bs-', linewidth=2, markersize=8, label='GPR (inverse)')
# O(N^3) 参考線
N_ref = np.array(Ns, dtype=float)
cubic_ref = times_krr[2] * (N_ref / N_ref[2])**3
ax.loglog(N_ref, cubic_ref, 'k--', alpha=0.5, label='$O(N^3)$ reference')
ax.set_xlabel('Number of training samples $N$', fontsize=12)
ax.set_ylabel('Training time (seconds)', fontsize=12)
ax.set_title('Computational Scaling of KRR and GPR', fontsize=14)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3, which='both')
plt.tight_layout()
plt.savefig('krr_computational_scaling.png', dpi=150, bbox_inches='tight')
plt.show()
計算時間は $N$ の 3 乗にほぼ比例して増加し、$N = 2000$ 程度ではまだ実用的ですが、$N = 10000$ を超えると近似手法が必要になることが確認できます。KRR と GPR は同じ $O(N^3)$ スケーリングを持ちますが、GPR は $\bm{K}^{-1}$ を保持するため、予測時の分散計算で追加のコストがかかります。
まとめ
本記事では、カーネルリッジ回帰と RKHS の理論を詳しく解説しました。
- リッジ回帰の双対形式: Woodbury の恒等式により、$d \times d$ の逆行列を $N \times N$ の逆行列に変換でき、高次元特徴空間で有利になります
- リプレゼンター定理: RKHS での正則化付き経験リスク最小化の最適解が $f^* = \sum_i \alpha_i^* k(\bm{x}_i, \cdot)$ の形で表されることを証明しました。鍵は、$\mathcal{V}^{\perp}$ 成分がデータ点での値に寄与せず、正則化項を増加させるだけという点です
- カーネルリッジ回帰の閉形式解: $\bm{\alpha}^* = (\bm{K} + \lambda\bm{I})^{-1}\bm{y}$ で、カーネル関数のみで計算が完結します
- ガウス過程回帰との等価性: $\lambda = \sigma_n^2$ のとき予測関数が一致します。RKHS ノルムによる正則化がガウス過程の事前分布に対応しますが、GPR は不確実性の定量化も可能です
- 計算量: 双対形式は $O(N^3)$ で、大規模データにはナイストレム近似やランダム特徴などの近似手法が必要です
次のステップとして、以下の記事も参考にしてください。