ガウス過程回帰メモ(2)

この記事はアドベントカレンダー ほぼ厚木の民 の11日目記事です。 昨年の同じ時期に書いて放置していたエントリの続き(といっても部分的な走り書きのメモ)を書きます。

前回のおさらい

  • y=wTϕ(x)y = \boldsymbol{w}^{\mathsf{T}}\boldsymbol{\phi}(x) のパラメータ w\boldsymbol{w} を求めたいが、次元が指数的に増える
  • wN(0,λ2I)\boldsymbol{w} \sim \mathcal{N}(\boldsymbol{0}, \lambda^2 \boldsymbol{\mathrm{I}}) として期待値をとりパラメータを積分消去
  • yN(0,λ2ΦΦT)\boldsymbol{y} \sim \mathcal{N}\left( \boldsymbol{0}, \lambda^2 \Phi \Phi^{\mathrm{T}} \right)

ガウス過程回帰モデル

N個の観測値

D={(x1,y1),(x2,y2),,(xN,yN)}\mathcal{D} = \lbrace \left(\boldsymbol{x}_1, y_1\right), \left(\boldsymbol{x}_2, y_2\right), \cdots, \left(\boldsymbol{x}_N, y_N\right) \rbrace

について

y=f(x)y = f(\boldsymbol{x})

の関係があり、この ff がガウス過程 GP(0,k(x,x))GP(\boldsymbol{0}, k(\boldsymbol{x}, \boldsymbol{x}^{\prime})) に従うものとする。

このとき、 D\mathcal{D} にない新しい観測値 x\boldsymbol{x}^* が与えられたとき、どのように yy^*を予測するかを考える。 x1,x2,xN,x\boldsymbol{x}_1, \boldsymbol{x}_2, \cdots \boldsymbol{x}_N, \boldsymbol{x^*} から計算したカーネル行列を KK^{\prime}(y1,y2,,yN,y)T(y_1, y_2, \cdots, y_N, y^*)^{\mathsf{T}}y\boldsymbol{y}^{\prime} とおく。y\boldsymbol{y}^{\prime} はガウス分布 N(0,K)\mathcal{N}(\boldsymbol{0}, K^{\prime}) に従う

( yy)N(0,( KkkTk)K)\begin{aligned} \displaystyle \left( \begin{array}{c}   \boldsymbol{y}\\ y^* \end{array} \right) \sim \mathcal{N} \left( \boldsymbol{0}, \underbrace { \left( \begin{array}{cc}   K & \boldsymbol{k}_*\\ \boldsymbol{k}_*^{\mathsf{T}} & k_{**} \end{array} \right) }_{K^{\prime}} \right) \end{aligned}

ので、yy^* の分布を求めることができる。

p(yx,D)=N(kTK1y,kkTK1k)p(y^* | \boldsymbol{x}^*, \mathcal{D}) = \mathcal{N}\left(\boldsymbol{k}_*^{\mathsf{T}}K^{-1}\boldsymbol{y}, k_{**}-\boldsymbol{k}_{*}^{\mathsf{T}}K^{-1}\boldsymbol{k}_{*} \right)

x\boldsymbol{x}^* が複数与えられた場合(XX^*)も同様で、

p(yX,D)=N(kTK1y,kkTK1k)p(\boldsymbol{y}^* | X^*, \mathcal{D}) = \mathcal{N}\left(\boldsymbol{k}_*^{\mathsf{T}}K^{-1}\boldsymbol{y}, \boldsymbol{k}_{**}-\boldsymbol{k}_{*}^{\mathsf{T}}K^{-1}\boldsymbol{k}_{*} \right)

となる。

カーネルのハイパーパラメータ推定

例えばRBFカーネル

k(x,x)=θ1exp((xx)2θ2)+θ3δ(x,x)k(\boldsymbol{x}, \boldsymbol{x}^{\prime}) = \theta_1 \exp \left(-\frac{|(\boldsymbol{x}-\boldsymbol{x}^{\prime})|^2}{\theta_2}\right) + \theta_3 \delta(\boldsymbol{x},\boldsymbol{x}^{\prime})

のパラメータ θ=(θ1,θ2,θ3)\boldsymbol{\theta}=(\theta_1, \theta_2, \theta_3) についていい感じに推定したい。このとき θ\boldsymbol{\theta} に依存するカーネル行列を KθK_{\boldsymbol{\theta}} とおく。学習データの確率は

p(yX,θ)=N(y0,Kθ)=1(2π)N21Kθ12exp(12yTKy1y)\begin{aligned} p(\boldsymbol{y} | X, \boldsymbol{\theta}) & = \mathcal{N}\left(\boldsymbol{y} | \boldsymbol{0}, K_{\boldsymbol{\theta}}\right)\\ & = \frac{1}{(2\pi)^\frac{N}{2}} \frac{1}{|K_{\boldsymbol{\theta}}|^{\frac{1}{2}}} \exp \left(-\frac{1}{2}\boldsymbol{y}^{\mathsf{T}}K_{\boldsymbol{y}}^{-1}\boldsymbol{y}\right) \end{aligned}

対数尤度 LL を最大化することを考える。

L=logp(yX,θ)logKθyTKθ1y+Const.L = \log p(\boldsymbol{y} | X, \boldsymbol{\theta}) \propto - \log |K_{\boldsymbol{\theta}}| -\boldsymbol{y}^{\mathsf{T}}K_{\boldsymbol{\theta}}^{-1}\boldsymbol{y} + \text{Const.}

θ\boldsymbol{\theta} のある要素 θ\theta について

Lθ=tr(Kθ1Kθθ)+(Kθ1y)TKθθ(Kθ1y)\frac{\partial L}{\partial \theta} = \mathrm{tr} \left(K_{\boldsymbol{\theta}}^{-1}\frac{\partial K_{\boldsymbol{\theta}}}{\partial\theta}\right) + (K_{\boldsymbol{\theta}}^{-1}\boldsymbol{y})^{\mathsf{T}} \frac{\partial K_{\boldsymbol{\theta}}}{\partial\theta} (K_{\boldsymbol{\theta}}^{-1}\boldsymbol{y})

のように微分することができ、最急降下法やL-BFGSなどを使って最適化することができる。


Written by@satopirka
A software engineer.

GitHubTwitter

© 2021 Minato Sato. All Rights Reserved