Gauss過程回帰モデルとBayes線形回帰の等価性についてのメモ書きです。
Gauss過程回帰
Gauss過程回帰 (Gaussian process regression) モデルは任意個数の関数値が多変量Gauss分布に従うと仮定する回帰モデル。平均ベクトルの各成分は平均関数 m、共分散行列の各成分はカーネル関数 k によって定まる。
t 番目の入力変数が \bm x^{(t)} \in \mathbb R^d、観測された応答が y^{(t)} \in \mathbb R とすると、n+1 個の y は以下の分布に従うとモデル化される。
\begin{aligned}
\begin{bmatrix}
y^{(1)} \\
\vdots \\
y^{(n)} \\
y^{(n+1)} \\
\end{bmatrix}
\sim
\mathcal N \!
\left(
\begin{bmatrix}
m(\bm x^{(1)}) \\
\vdots \\
m(\bm x^{(n)}) \\
m(\bm x^{(n+1)}) \\
\end{bmatrix},
\begin{bmatrix}
k(\bm x^{(1)}, \bm x^{(1)}) + \sigma^2 & \cdots & k(\bm x^{(1)}, \bm x^{(n)}) & k(\bm x^{(1)}, \bm x^{(n+1)}) \\
\vdots & \ddots & \vdots & \vdots \\
k(\bm x^{(n)}, \bm x^{(1)}) & \cdots & k(\bm x^{(n)}, \bm x^{(n)}) + \sigma^2 & k(\bm x^{(n)}, \bm x^{(n+1)}) \\
k(\bm x^{(n+1)}, \bm x^{(1)}) & \cdots & k(\bm x^{(n+1)}, \bm x^{(n)}) & k(\bm x^{(n+1)}, \bm x^{(n+1)}) + \sigma^2 \\
\end{bmatrix}
\right).
\end{aligned}
あるいは \bm y_{n} = [y^{(t)}]_{t} \in \mathbb R^{n}, \bm m_n = [m(\bm x^{(t)})]_{t} \in \mathbb R^{n}, \bm K_{n} = [k(\bm x^{(i)}, \bm x^{(j)})]_{i,j} \in \mathbb R^{n \times n}, \bm k_{n} = [k(\bm x^{(i)}, \bm x^{(n+1)})]_{i} \in \mathbb R^{n} とおけば
\begin{aligned}
\begin{bmatrix}
\bm y \\
y^{(n+1)}
\end{bmatrix}
\sim
\mathcal N \!
\left(
\begin{bmatrix}
\bm m_n \\
m(\bm x^{(n+1)})
\end{bmatrix},
\begin{bmatrix}
\bm K_{n} + \sigma^2 \bm I_{n} & \bm k_{n} \\
\bm k_{n}^\top & k(\bm x^{(n+1)}, \bm x^{(n+1)}) + \sigma^2
\end{bmatrix}
\right).
\end{aligned}
ここで共分散行列の対角成分に \sigma^2 が加わっているのは、観測値にノイズが含まれると仮定していることを意味する。
\begin{aligned}
y^{(t)} = f(\bm x^{(t)}) + \varepsilon^{(t)}, \quad \varepsilon^{(t)} \sim \mathcal N(0, \sigma^2).
\end{aligned}
\bm X_n = [x^{(t)}_i]_{t,i} \in \mathbb R^{n \times d}, \bm y_n が与えられた条件下の y^{(n+1)} の条件付き分布 (予測分布) は
\begin{aligned}
y^{(n+1)} \mid \bm x^{(n+1)}, \bm X_n, \bm y_n
&\sim \mathcal N \!
\left(
\mu_n, \sigma^2_n
\right),
\end{aligned}
\\
\left\{
\begin{aligned}
\mu_n &= m(\bm x^{(n+1)}) + \bm k_n^\top (\bm K_n + \sigma^2 \bm I_n)^{-1} (\bm y_n - \bm m_n), \\
\sigma^2_n &= k(\bm x^{(n+1)}, \bm x^{(n+1)}) + \sigma^2 - \bm k_n^\top (\bm K_n + \sigma^2 \bm I_n)^{-1} \bm k_n.
\end{aligned}
\right.
多くの場合は、\bm y_n が零平均になるように平行移動されているとして、m = 0 として扱い、
\begin{aligned}
y^{(n+1)} \mid \bm x^{(n+1)}, \bm X_n, \bm y_n
&\sim \mathcal N \!
\left(
\mu_n, \sigma^2_n
\right),
\end{aligned}
\\
\left\{
\begin{aligned}
\mu_n &= \bm k_n^\top (\bm K_n + \sigma^2 \bm I_n)^{-1} \bm y_n, \\
\sigma^2_n &= k(\bm x^{(n+1)}, \bm x^{(n+1)}) + \sigma^2 - \bm k_n^\top (\bm K_n + \sigma^2 \bm I_n)^{-1} \bm k_n
\end{aligned}
\right.
とする。さらに、場合によっては \bm y_n が与えられた条件下の “ノイズを含む” 観測値 y^{(n+1)} = f(\bm x^{(n+1)}) + \varepsilon の予測分布ではなく、“ノイズが生じる前の” 関数値 f^{(n+1)} = f(\bm x^{(n+1)}) そのものの予測分布を
\begin{gathered}
\begin{aligned}
f^{(n+1)} \mid \bm x^{(n+1)}, \bm X_n, \bm y_n
&\sim \mathcal N \!
\left(
\mu_n, \sigma^2_n
\right),
\end{aligned}
\\
\left\{
\begin{aligned}
\mu_n &= \bm k_n^\top (\bm K_n + \sigma^2 \bm I_n)^{-1} \bm y_n, \\
\sigma^2_n &= k(\bm x^{(n+1)}, \bm x^{(n+1)}) + \bm k_n^\top (\bm K_n + \sigma^2 \bm I_n)^{-1} \bm k_n
\end{aligned}
\right.
\end{gathered}
と求めることもある。
線形回帰との関係
Gauss過程回帰モデルは、一般には y^{(n+1)} や f^{(n+1)} を何らかの閉形式で書かれた関数 \hat f(\bm x^{(n+1)}; \bm \theta) を用いて表すことができるとは限らない。しかしカーネル関数がある条件を満たす場合には、y^{(n+1)} を線形回帰モデルで近似することができる。
平均関数が零の場合の y に関する予測分布
\begin{gathered}
\begin{aligned}
y^{(n+1)} \mid \bm x^{(n+1)}, \bm X_n, \bm y_n
&\sim \mathcal N \!
\left(
\mu_n, \sigma^2_n
\right),
\end{aligned}
\\
\begin{aligned}
\mu_n &= \bm k_n^\top (\bm K_n + \sigma^2 \bm I_n)^{-1} \bm y_n, \\
\sigma^2_n &= k(\bm x^{(n+1)}, \bm x^{(n+1)}) + \sigma^2 - \bm k_n^\top (\bm K_n + \sigma^2 \bm I_n)^{-1} \bm k_n
\end{aligned}
\end{gathered}
において
\begin{aligned}
\bm K_n
&=
\sigma^2
\bm X_n \bm \Sigma_0^{-1} \bm X_n^\top, \\
\bm k_n
&=
\sigma^2
\bm X_n \bm \Sigma_0^{-1} \bm x^{(n+1)}, \\
k(\bm x^{(n+1)}, \bm x^{(n+1)})
&=
\sigma^2
\bm x^{(n+1) \top} \bm \Sigma_0^{-1} \bm x^{(n+1)}
\end{aligned}
とすれば、Woodburyの恒等式を使って
\begin{aligned}
\mu_n
&=
\sigma^2
\bm x^{(n+1) \top} \bm \Sigma_0^{-1} \bm X_n^\top
\left(
\sigma^2
\bm X_n
\bm \Sigma_0^{-1}
\bm X_n^\top
+
\sigma^2 \bm I_n
\right)^{-1} \bm y_n
\\
&=
\bm x^{(n+1) \top} \bm \Sigma_0^{-1} \bm X_n^\top
\left(
\bm X_n
\bm \Sigma_0^{-1}
\bm X_n^\top
+
\bm I_n
\right)^{-1} \bm y_n
\\
&=
\bm x^{(n+1) \top} \bm \Sigma_0^{-1} \bm X_n^\top
\left(
\bm I_n
-
\left(
\bm I_n
+
\bm X_n \bm \Sigma_0^{-1} \bm X_n^\top
\right)^{-1}
\bm X_n \bm \Sigma_0^{-1} \bm X_n^\top
\right) \bm y_n
\\
&=
\bm x^{(n+1) \top}
\left(
\bm \Sigma_0^{-1}
-
\bm \Sigma_0^{-1} \bm X_n^\top
\left(
\bm I_n
+
\bm X_n \bm \Sigma_0^{-1} \bm X_n^\top
\right)^{-1}
\bm X_n
\bm \Sigma_0^{-1}
\right)
\bm X_n^\top \bm y_n
\\
&=
\bm x^{(n+1) \top} \underbrace{
\left(
\bm X_n^\top \bm X_n + \bm \Sigma_0
\right)^{-1}
\bm X_n^\top
\bm y_n
}_{\text{post. mean of } \bm w}.
\end{aligned}
同様に
\begin{aligned}
\sigma^2_n
&=
\sigma^2
\bm x^{(n+1) \top} \bm \Sigma_0^{-1} \bm x^{(n+1)}
+
\sigma^2
\\&\qquad
-
\sigma^2
\bm x^{(n+1) \top} \bm \Sigma_0^{-1} \bm X_n^\top
\left(
\sigma^2
\bm X_n \bm \Sigma_0^{-1} \bm X_n^\top
+
\sigma^2 \bm I_n
\right)^{-1}
\sigma^2
\bm X_n \bm \Sigma_0^{-1} \bm x^{(n+1)}
\\
&=
\sigma^2
\bm x^{(n+1) \top}
\left(
\bm \Sigma_0^{-1}
-
\bm \Sigma_0^{-1} \bm X_n^\top
\left(
\bm I_n
+
\bm X_n \bm \Sigma_0^{-1} \bm X_n^\top
\right)^{-1}
\bm X_n
\bm \Sigma_0^{-1}
\right)
\bm x^{(n+1)}
\\&\qquad
+
\sigma^2
\\
&=
\bm x^{(n+1) \top}
\underbrace{
\sigma^2
\left( \bm X_n^\top \bm X_n + \bm \Sigma_0^{-1} \right)^{-1}
}_{\text{post. cov. of } \bm w}
\bm x^{(n+1)} + \sigma^2.
\end{aligned}
これは線形回帰モデル
\begin{aligned}
y &= \bm x^\top \bm w + \varepsilon, \\
\varepsilon &\sim \mathcal N(0, \sigma^2)
\end{aligned}
において、\bm w の事前分布として
\begin{aligned}
\bm w &\sim \mathcal N_d(\bm 0, \sigma^2 \bm \Sigma_0)
\end{aligned}
を仮定したときの事後予測分布に一致する。実際、\bm w の事後分布は正規分布となって
\begin{aligned}
\bm w \mid \bm X_n, \bm y_n
&\sim
\mathcal N_d
\left(
\bm A_n^{-1} \bm X_n^\top \bm y_n,
\sigma^2 \bm A_n^{-1}
\right),
\end{aligned}
\\
\begin{aligned}
\bm A_n &= \bm X_n^\top \bm X_n + \bm \Sigma_0^{-1}
\end{aligned}
となるから、このときの y^{(n+1)} の予測分布は正規分布となる。平均と分散は
\begin{aligned}
\mathbb E_{y^{(n+1)} \mid \bm X_n, \bm y_n} [y^{(n+1)}]
&=
\mathbb E_{\bm w \mid \bm X_n, \bm y_n} [\bm x^{(n+1) \top} \bm w]
\\
&=
\bm x^{(n+1) \top} \mathbb E_{\bm w \mid \bm X_n, \bm y_n} [\bm w]
\\
&=
\bm x^{(n+1) \top} \bm A_n^{-1} \bm X_n^\top \bm y_n,
\\
\mathbb V_{y^{(n+1)} \mid \bm X_n, \bm y_n} [y^{(n+1)}]
&=
\mathbb V_{\bm w \mid \bm X_n, \bm y_n} [\bm x^{(n+1) \top} \bm w]
+ \sigma^2
\\
&=
\bm x^{(n+1) \top} \mathrm{Cov}_{\bm w \mid \bm X_n, \bm y_n} [\bm w] \bm x^{(n+1)}
+ \sigma^2
\\
&=
\bm x^{(n+1) \top} \sigma^2 \bm A_n^{-1} \bm x^{(n+1)} + \sigma^2.
\end{aligned}
Weight space view
このように、何らかの正定値行列 \bm \Sigma_0 をベクトルで挟むようにしてカーネル関数を
\begin{aligned}
k(\bm x, \bm x') = \sigma^2 \bm x^\top \bm \Sigma_0^{-1} \bm x'
\end{aligned}
と定めると、Gauss過程回帰モデルの予測分布は線形回帰モデル
\begin{aligned}
y &= \bm x^\top \bm w + \varepsilon, \\
\varepsilon &\sim \mathcal N(0, \sigma^2), \\
\bm w &\sim \mathcal N_d(\bm 0, \sigma^2 \bm \Sigma_0)
\end{aligned}
の予測分布に一致する。
\bm x \in \mathbb R^d を特徴空間 \bm \phi \in \mathbb R^p にマッピングして
\begin{aligned}
k(\bm x, \bm x') = \sigma^2 \bm \phi(\bm x)^\top \bm \Sigma_0^{-1} \bm \phi(\bm x')
\end{aligned}
とする場合は、線形回帰モデルとしては
\begin{aligned}
y &= \bm \phi(\bm x)^\top \bm \theta + \varepsilon, \\
\varepsilon &\sim \mathcal N(0, \sigma^2), \\
\bm \theta &\sim \mathcal N_p(\bm 0, \sigma^2 \bm \Sigma_0)
\end{aligned}
を考えていると思えば良い。
\bm \phi が無限次元だと考えて、
\begin{gathered}
k(\bm x, \bm x')
= \sigma^2 \sigma_0^2 \bm \phi(\bm x)^\top \bm \phi(\bm x')
= \sigma^2 \sigma_0^2 r(\bm x, \bm x'), \\
r(\bm x, \bm x') = \sum_{i=1}^\infty \phi_i(\bm x) \phi_i(\bm x')
\end{gathered}
というようなカーネルを考えることが多い。
しばしば weight space view、日本語にすると重み空間の視点とか重み空間表示とかいう用語が出てくることがあるが、これは以上のようにカーネル関数を何らかの特徴量の内積と捉えてGauss過程回帰モデルを線形回帰モデルとしてみなす見方のことを指している。
Discussion