6/27 追記: typo, p(y∣f) p(\mathbf{y} | \mathbf{f}) の誤字を修正, qFITC(f∗∣y) q_{\text{FITC}}(\mathbf{f}_* | \mathbf{y}) の二番目の等号を修正 (σ−2 \sigma^{-2} を削除)

「ガウス過程と機械学習」を読んでいるが、5.2補助変数法のところで、どの部分で近似が行われているのかよく分からなくなってしまった。

そのため、今回は原論文であるQuinonero Candela, J. and Rasmussen, CE.の “A Unifying View of Sparse Approximate Gaussian Process Regression” を読んでスパース近似についてまとめて見ようと思う。ゴールは、The Fully Independent Training Conditional (FITC) の理解である。

X=(x1,…,xN) \mathbf{X}=(\mathbf{x}_1, \ldots, \mathbf{x}_N) を学習データ、 y=(y1,…,yN)⊤ \mathbf{y}=(y_1, \ldots, y_N)^\top を観測値とする。学習データと観測値の関係は、ガウス過程から生成される関数 f f と誤差 ϵn \epsilon_n を用いて

yn=f(xn)+ϵn, y_n = f(\mathbf{x_n}) + \epsilon_n,

ϵn∼N(0,σ2)\epsilon_n \sim \mathcal{N}(0, \sigma^2)

と結びつく観測モデルを考える。f=(f1,…,fN)⊤,fn=f(xn) \mathbf{f}=(f_1, \ldots, f_N)^\top, f_n=f(\mathbf{x}_n) は学習データの出力値である。 上の観測モデルは、

p(y∣f)=N(f,σ2I) p(\mathbf{y} | \mathbf{f}) = \mathcal{N}(\mathbf{f}, \sigma^2\mathbf{I})

と書き直せる。

取り組みたい問題は、ガウス過程回帰に基づいた回帰モデル。

予測したい点を X∗=(x∗1,…,x∗M) \mathbf{X}_*=(\mathbf{x}_{*1}, \ldots, \mathbf{x}_{*M} ) , 出力値を f∗ \mathbf{f}_* , 観測値を y∗ \mathbf{y}_* とする。

f,f∗ \mathbf{f}, \mathbf{f}_* の条件付き同時確率分布はベイズルールから

p(f,f∗∣y)=p(f,f∗)p(f∣y)p(y), \begin{aligned} p(\mathbf{f}, \mathbf{f}_* | \mathbf{y}) = \frac{p(\mathbf{f}, \mathbf{f}_*)p(\mathbf{f} | \mathbf{y})}{p(\mathbf{y})}, \end{aligned}

同時事前確率は、ガウス過程の定義より、カーネル関数が定める共分散行列 K=(k(xi,xj))i,j \mathbf{K}=(k(\mathbf{x}_i, \mathbf{x}_j))_{i, j} を用いて

p(f,f∗)=N(0,(Kf,fK∗,fKf,∗K∗,∗)), \begin{aligned} p(\mathbf{f}, \mathbf{f}_*) = \mathcal{N}\left( \mathbf{0}, \left( \begin{matrix} \mathbf{K}_{\mathbf{f}, \mathbf{f}} & \mathbf{K}_{*, \mathbf{f}} \\ \mathbf{K}_{\mathbf{f}, *} & \mathbf{K}_{*, *} \end{matrix} \right) \right), \end{aligned}

となる。学習データ X \mathbf{X} が大きくなると Kf,f \mathbf{K}_{\mathbf{f}, \mathbf{f}} の計算量が大きくなるので、この部分の計算量を減らすような近似を行いたいわけである。

ガウス過程の予測分布 p(f∗∣y) p(\mathbf{f}_* | \mathbf{y} ) は、

p(f∗∣y)=∫p(f,f∗∣y)df=1p(y)∫p(y∣f)p(f,f∗)df=N(K∗,f K~f,f−1 y,K∗,∗−K∗,f K~f,f−1 Kf,∗), \begin{aligned} p(\mathbf{f}_* | \mathbf{y}) &= \int p (\mathbf{f}, \mathbf{f}_* | \mathbf{y}) d\mathbf{f} \\ &= \frac{1}{p(\mathbf{y})} \int p(\mathbf{y} | \mathbf{f}) p(\mathbf{f}, \mathbf{f}_*)d\mathbf{f} \\ &= \mathcal{N}(\mathbf{K}_{*, \mathbf{f}}\ \widetilde{\mathbf{K}}_{\mathbf{f}, \mathbf{f}}^{-1}\ \mathbf{y}, \mathbf{K}_{*, *} - \mathbf{K}_{*, \mathbf{f}}\ \widetilde{\mathbf{K}}_{\mathbf{f}, \mathbf{f}}^{-1}\ \mathbf{K}_{\mathbf{f}, *}), \end{aligned}

ここで K~f,f=Kf,f+σ2I \widetilde{\mathbf{K}}_{\mathbf{f}, \mathbf{f}} = \mathbf{K}_{\mathbf{f}, \mathbf{f}} + \sigma^2 \mathbf{I} である。 この式は「ガウス過程と機械学習」の公式3.8に対応し、こちらではノイズ項目 σ2I \sigma^2 \mathbf{I} もカーネルに含めているので、notationが少し違っている。

いくつかスパース近似にはバリエーションがあるが、補助入力点 u=(u1,…,um)⊤ \mathbf{u}=(\mathbf{u}_1, \ldots, \mathbf{u}_m)^\top を使って p(f∗∣f) p(\mathbf{f}_* | \mathbf{f}) を近似するという方法は共通している。まずは近似ではなく正確に成り立っている式を確認する。

p(f∗,f)=∫p(f∗,f,u)du=∫p(f∗,f∣u)p(u)du, \begin{aligned} p(\mathbf{f}_*, \mathbf{f}) &= \int p(\mathbf{f}_*, \mathbf{f}, \mathbf{u})d\mathbf{u} \\ &= \int p(\mathbf{f}_*, \mathbf{f} | \mathbf{u})p(\mathbf{u})d\mathbf{u}, \end{aligned} p(u)=N(0,Ku,u), p(\mathbf{u}) = \mathcal{N}(\mathbf{0}, \mathbf{K}_{\mathbf{u}, \mathbf{u}}), p(f∣u)=N(Kf,uKu,u−1u,Kf,f−Qf,f), p(\mathbf{f} | \mathbf{u}) = \mathcal{N}(\mathbf{K}_{\mathbf{f}, \mathbf{u}} \mathbf{K}_{\mathbf{u}, \mathbf{u}}^{-1} \mathbf{u}, \mathbf{K}_{\mathbf{f}, \mathbf{f}} - \mathbf{Q}_{\mathbf{f}, \mathbf{f}}), p(f∗∣u)=N(K∗,uKu,u−1u,K∗,∗−Q∗,∗) p(\mathbf{f}_* | \mathbf{u}) = \mathcal{N}(\mathbf{K}_{*, \mathbf{u}} \mathbf{K}_{\mathbf{u}, \mathbf{u}}^{-1} \mathbf{u}, \mathbf{K}_{*, *} - \mathbf{Q}_{*, *})

と表される。ここで、 Qa,b:=Ka,uKu,u−1Ku,b \mathbf{Q}_{\mathbf{a}, \mathbf{b}} := \mathbf{K}_{\mathbf{a}, \mathbf{u}} \mathbf{K}_{\mathbf{u}, \mathbf{u}}^{-1} \mathbf{K}_{\mathbf{u}, \mathbf{b}} である。

補助入力点 u \mathbf{u} を使った条件付き確率 q(f∗∣u),q(f∣u) q(\mathbf{f}_* | \mathbf{u}), q(\mathbf{f} | \mathbf{u}) を使って、

p(f∗,f)≃q(f∗,f)=∫q(f∗∣u)q(f∣u)p(u)du \begin{aligned} p(\mathbf{f}_*, \mathbf{f}) &\simeq q(\mathbf{f}_*, \mathbf{f}) \\ &= \int q(\mathbf{f}_* | \mathbf{u}) q(\mathbf{f} | \mathbf{u}) p(\mathbf{u}) d\mathbf{u} \end{aligned}

という近似を行う。この時の q q の作り方が近似の手法によって異なることになる。

The Subset of Data (SoD) Approximation

部分データ法は元々の入力点の中から代表となる点を抜き出し、抜き出した点を新たな入力点としてガウス回帰を行うやり方。細かい説明は省略。

The Subset of Regressors (SoR) Approximation

The Subset of Regressors Approximation (部分回帰法?, 以下SoR) は f,f∗ \mathbf{f}, \mathbf{f}_* の同時分布を

qSoR(f,f∗)=N(0,(Qf,fQ∗,fQf,∗Q∗,∗)) \begin{aligned} q_{\text{SoR}}(\mathbf{f}, \mathbf{f}_*) = \mathcal{N}\left( \mathbf{0}, \left( \begin{matrix} \mathbf{Q}_{\mathbf{f}, \mathbf{f}} & \mathbf{Q}_{*, \mathbf{f}} \\ \mathbf{Q}_{\mathbf{f}, *} & \mathbf{Q}_{*, *} \end{matrix} \right) \right) \end{aligned}

と近似するもので、後に出てくるDTC, FITCの基礎となる考え方である。 心持ちとしては、

qSoR(f∣u)=N(Kf,uKu,u−1u,0)≃N(Kf,uKu,u−1u,Kf,f−Qf,f)=p(f∣u), \begin{aligned} q_{\text{SoR}}(\mathbf{f} | \mathbf{u}) &= \mathcal{N}(\mathbf{K}_{\mathbf{f}, \mathbf{u}} \mathbf{K}_{\mathbf{u}, \mathbf{u}}^{-1} \mathbf{u}, \mathbf{0}) \\ &\simeq \mathcal{N}(\mathbf{K}_{\mathbf{f}, \mathbf{u}} \mathbf{K}_{\mathbf{u}, \mathbf{u}}^{-1} \mathbf{u}, \mathbf{K}_{\mathbf{f}, \mathbf{f}} - \mathbf{Q}_{\mathbf{f}, \mathbf{f}}) \\ &= p(\mathbf{f} | \mathbf{u}), \end{aligned} qSoR(f∗∣u)=N(K∗,uKu,u−1u,0)≃N(K∗,uKu,u−1u,K∗,∗−Q∗,∗)=p(f∗∣u) \begin{aligned} q_{\text{SoR}}(\mathbf{f}_* | \mathbf{u}) &= \mathcal{N}(\mathbf{K}_{*, \mathbf{u}} \mathbf{K}_{\mathbf{u}, \mathbf{u}}^{-1} \mathbf{u}, \mathbf{0}) \\ &\simeq \mathcal{N}(\mathbf{K}_{*, \mathbf{u}} \mathbf{K}_{\mathbf{u}, \mathbf{u}}^{-1} \mathbf{u}, \mathbf{K}_{*, *} - \mathbf{Q}_{*, *}) \\ &= p(\mathbf{f}_* | \mathbf{u}) \end{aligned}

と、p(f∣u),p(f∗∣u) p(\mathbf{f} | \mathbf{u}), p(\mathbf{f}_* | \mathbf{u}) の分散を 0 0 と近似して、f,f∗ \mathbf{f}, \mathbf{f}_* と u \mathbf{u} の関係がdeterministicと仮定するものである。(deterministicの場合はもはや f∣u \mathbf{f} | \mathbf{u} は正規分布ではないので、この書き方は微妙な感じかもしれないが、論文のnotationに従った)

仮定から qSoR(f,f∗) q_{\text{SoR}}(\mathbf{f}, \mathbf{f}_* ) の導出は、 同時分布が正規分布になることと、多変数正規分布の線形変換の性質から、

f=Kf,uKu,u−1u∼N(0,(Kf,uKu,u−1)Ku,u(Kf,uKu,u−1)⊤)∼N(0,Kf,uKu,u−1Ku,f)∼N(0,Qf,f), \begin{aligned} \mathbf{f} &= \mathbf{K}_{\mathbf{f}, \mathbf{u}} \mathbf{K}_{\mathbf{u}, \mathbf{u}}^{-1} \mathbf{u} \\ & \sim \mathcal{N}(\mathbf{0}, (\mathbf{K}_{\mathbf{f}, \mathbf{u}} \mathbf{K}_{\mathbf{u}, \mathbf{u}}^{-1}) \mathbf{K}_{\mathbf{u}, \mathbf{u}} (\mathbf{K}_{\mathbf{f}, \mathbf{u}} \mathbf{K}_{\mathbf{u}, \mathbf{u}}^{-1})^\top) \\ & \sim \mathcal{N}(\mathbf{0}, \mathbf{K}_{\mathbf{f}, \mathbf{u}} \mathbf{K}_{\mathbf{u}, \mathbf{u}}^{-1} \mathbf{K}_{\mathbf{u}, \mathbf{f}}) \\ & \sim \mathcal{N}(\mathbf{0}, \mathbf{Q}_{\mathbf{f}, \mathbf{f}}), \end{aligned}

であり、同様に f∗∼N(0,Q∗,∗) \mathbf{f}_* \sim \mathcal{N}(\mathbf{0}, \mathbf{Q}_{*, *}) を得られることと、相互共分散行列の線型性から、

cov(f,f∗)=cov(Kf,uKu,u−1u,K∗,uKu,u−1u)=Kf,uKu,u−1cov(u,u)(K∗,uKu,u−1)⊤=Kf,uKu,u−1Ku,u−1Ku,∗=Qf,∗, \begin{aligned} \text{cov}(\mathbf{f}, \mathbf{f}_*) &= \text{cov}(\mathbf{K}_{\mathbf{f}, \mathbf{u}} \mathbf{K}_{\mathbf{u}, \mathbf{u}}^{-1} \mathbf{u}, \mathbf{K}_{*, \mathbf{u}} \mathbf{K}_{\mathbf{u}, \mathbf{u}}^{-1} \mathbf{u} ) \\ & = \mathbf{K}_{\mathbf{f}, \mathbf{u}} \mathbf{K}_{\mathbf{u}, \mathbf{u}}^{-1} \text{cov}(\mathbf{u}, \mathbf{u}) (\mathbf{K}_{*, \mathbf{u}} \mathbf{K}_{\mathbf{u}, \mathbf{u}}^{-1})^\top \\ & = \mathbf{K}_{\mathbf{f}, \mathbf{u}} \mathbf{K}_{\mathbf{u}, \mathbf{u}}^{-1} \mathbf{K}_{\mathbf{u}, \mathbf{u}}^{-1} \mathbf{K}_{\mathbf{u}, *} \\ & = \mathbf{Q}_{\mathbf{f}, *}, \end{aligned}

cov(f∗,f)=Q∗,f \text{cov}(\mathbf{f}_*, \mathbf{f}) = \mathbf{Q}_{*, \mathbf{f}} であることから従う。

予測分布は、

qSoR(f∗∣y)=N(Q∗,f(Qf,f+σ2I)−1y,Q∗,∗−Q∗,f(Qf,f+σ2I)−1Qf,∗) \begin{aligned} q_{\text{SoR}}(\mathbf{f}_* | \mathbf{y}) & = \mathcal{N}(\mathbf{Q}_{*, \mathbf{f}}(\mathbf{Q}_{\mathbf{f},\mathbf{f}} + \sigma^2 \mathbf{I})^{-1} \mathbf{y}, \mathbf{Q}_{*, *} - \mathbf{Q}_{*, \mathbf{f}}(\mathbf{Q}_{\mathbf{f},\mathbf{f}} + \sigma^2 \mathbf{I})^{-1} \mathbf{Q}_{\mathbf{f}, *}) \end{aligned}

で与えられる。これが Σ=(σ−2Ku,fKf,u+Ku,u)−1 \Sigma = (\sigma^{-2} \mathbf{K}_{\mathbf{u}, \mathbf{f}} \mathbf{K}_{\mathbf{f}, \mathbf{u}} + \mathbf{K}_{\mathbf{u}, \mathbf{u}})^{-1} を用いて N(σ−2K∗,uΣKu,fy,K∗,uΣKu,∗ \mathcal{N} (\sigma^{-2} \mathbf{K}_{*, \mathbf{u}} \Sigma \mathbf{K}_{\mathbf{u}, \mathbf{f}} \mathbf{y}, \mathbf{K}_{*, \mathbf{u}} \Sigma \mathbf{K}_{\mathbf{u}, *} ) と等しくなることを示そう。 初見でなぜこうなるか不明だった……

まず、

σ2Σ−1Ku,u−1Ku,f=(Ku,fKf,u+σ2Ku,u)Ku,u−1Ku,f=Ku,f(Kf,uKu,u−1Ku,f+σ2I)=Ku,f(Qf,f+σ2I) \begin{aligned} \sigma^2 \Sigma^{-1} \mathbf{K}_{\mathbf{u}, \mathbf{u}}^{-1} \mathbf{K}_{\mathbf{u}, \mathbf{f}} & = (\mathbf{K}_{\mathbf{u}, \mathbf{f}} \mathbf{K}_{\mathbf{f}, \mathbf{u}} + \sigma^2 \mathbf{K}_{\mathbf{u}, \mathbf{u}}) \mathbf{K}_{\mathbf{u}, \mathbf{u}}^{-1} \mathbf{K}_{\mathbf{u}, \mathbf{f}} \\ & = \mathbf{K}_{\mathbf{u}, \mathbf{f}} (\mathbf{K}_{\mathbf{f}, \mathbf{u}} \mathbf{K}_{\mathbf{u}, \mathbf{u}}^{-1} \mathbf{K}_{\mathbf{u}, \mathbf{f}} + \sigma^2 \mathbf{I}) \\ & = \mathbf{K}_{\mathbf{u}, \mathbf{f}} (\mathbf{Q}_{\mathbf{f}, \mathbf{f}} + \sigma^2 \mathbf{I}) \end{aligned}

であるから、

Q∗,f(Qf,f+σ2I)−1=σ−2K∗,uΣ(σ2Σ−1Ku,u−1Ku,f)(Qf,f+σ2I)−1=σ−2K∗,uΣKu,f \begin{aligned} \mathbf{Q}_{*, \mathbf{f}}(\mathbf{Q}_{\mathbf{f},\mathbf{f}} + \sigma^2 \mathbf{I})^{-1} & = \sigma^{-2} \mathbf{K}_{*, \mathbf{u}} \Sigma (\sigma^2 \Sigma^{-1} \mathbf{K}_{\mathbf{u}, \mathbf{u}}^{-1} \mathbf{K}_{\mathbf{u}, \mathbf{f}}) (\mathbf{Q}_{\mathbf{f}, \mathbf{f}} + \sigma^2 \mathbf{I})^{-1} \\ & = \sigma^{-2} \mathbf{K}_{*, \mathbf{u}} \Sigma \mathbf{K}_{\mathbf{u}, \mathbf{f}} \end{aligned}

を得る。分散の方の等号も、この式を使えば最終的に Σ(Ku,u+σ−2Ku,fKf,u)=I \Sigma (\mathbf{K}_{\mathbf{u}, \mathbf{u}} + \sigma^{-2}\mathbf{K}_{\mathbf{u}, \mathbf{f}} \mathbf{K}_{\mathbf{f}, \mathbf{u}}) = \mathbf{I} に帰着させて示すことができる。

The Deterministic Training Conditional (DTC) Approximation

DTCはSoRと似ていて、

qDTC(f∣u)=N(Kf,uKu,u−1u,0) q_{\text{DTC}}(\mathbf{f} | \mathbf{u}) = \mathcal{N}(\mathbf{K}_{\mathbf{f}, \mathbf{u}} \mathbf{K}_{\mathbf{u}, \mathbf{u}}^{-1} \mathbf{u}, \mathbf{0})

と f \mathbf{f} が deterministic に u \mathbf{u} によって定まると仮定する。一方、出力値の方は

qDTC(f∗∣u)=p(f∗∣u)=N(K∗,uKu,u−1u,K∗,∗−Q∗,∗) \begin{aligned} q_{\text{DTC}}(\mathbf{f}_* | \mathbf{u}) &= p(\mathbf{f}_* | \mathbf{u}) \\ &= \mathcal{N}(\mathbf{K}_{*, \mathbf{u}} \mathbf{K}_{\mathbf{u}, \mathbf{u}}^{-1} \mathbf{u}, \mathbf{K}_{*, *} - \mathbf{Q}_{*, *}) \end{aligned}

と近似を行わない。

cov(f,f∗)=K∗,uKu,u−1cov(u,f∗)=K∗,uKu,u−1Ku,∗=Qf,∗ \begin{aligned} \text{cov}(\mathbf{f}, \mathbf{f}_*) &= \mathbf{K}_{*, \mathbf{u}} \mathbf{K}_{\mathbf{u}, \mathbf{u}}^{-1} \text{cov}(\mathbf{u}, \mathbf{f}_*) \\ &= \mathbf{K}_{*, \mathbf{u}} \mathbf{K}_{\mathbf{u}, \mathbf{u}}^{-1} \mathbf{K}_{\mathbf{u}, *} \\ &= \mathbf{Q}_{\mathbf{f}, *} \end{aligned}

などから、

qDTC(f,f∗)=N(0,(Qf,fQ∗,fQf,∗K∗,∗)) \begin{aligned} q_{\text{DTC}}(\mathbf{f}, \mathbf{f}_*) = \mathcal{N}\left( \mathbf{0}, \left( \begin{matrix} \mathbf{Q}_{\mathbf{f}, \mathbf{f}} & \mathbf{Q}_{*, \mathbf{f}} \\ \mathbf{Q}_{\mathbf{f}, *} & \mathbf{K}_{*, *} \end{matrix} \right) \right) \end{aligned}

となる。SoRの時とほとんど同じ。入力値の近似で計算量が十分削減できる場合は、出力値のカーネルの部分で正確な値を使いたいということだろう。

予測分布は、

qDTC(f∗∣y)=N(Q∗,f(Qf,f+σ2I)−1y,K∗,∗−Q∗,f(Qf,f+σ2I)−1Qf,∗)=N(σ−2K∗,uΣKu,fy,K∗,∗−Q∗,∗+K∗,uΣKu,∗) \begin{aligned} q_{\text{DTC}}(\mathbf{f}_* | \mathbf{y}) & = \mathcal{N}(\mathbf{Q}_{*, \mathbf{f}}(\mathbf{Q}_{\mathbf{f},\mathbf{f}} + \sigma^2 \mathbf{I})^{-1} \mathbf{y}, \mathbf{K}_{*, *} - \mathbf{Q}_{*, \mathbf{f}}(\mathbf{Q}_{\mathbf{f},\mathbf{f}} + \sigma^2 \mathbf{I})^{-1} \mathbf{Q}_{\mathbf{f}, *}) \\ & = \mathcal{N} (\sigma^{-2} \mathbf{K}_{*, \mathbf{u}} \Sigma \mathbf{K}_{\mathbf{u}, \mathbf{f}} \mathbf{y}, \mathbf{K}_{*, *} - \mathbf{Q}_{*, *} + \mathbf{K}_{*, \mathbf{u}} \Sigma \mathbf{K}_{\mathbf{u}, *}) \end{aligned}

となる。これはSoRの時の結果からすぐに従う。

The Fully Independent Training Conditional (FITC) Approximation

いよいよ、当初の目的だったFITCまでたどり着いた。FITCは、DTCとほぼ同じだが、 DTCでは切り捨てていた q(f∣u) q(\mathbf{f} | \mathbf{u}) の分散を考慮する。

ただ、分散共分散行列をフルで計算したくない (フルで計算すると近似を行わない通常のガウス回帰である) ので、 出力値の間の相関を無視して、対角線の部分だけを計算する。つまり、

qFITC(f∣u)=N(Kf,uKu,u−1u,diag(Kf,f−Qf,f)) q_{\text{FITC}}(\mathbf{f} | \mathbf{u}) = \mathcal{N}(\mathbf{K}_{\mathbf{f}, \mathbf{u}} \mathbf{K}_{\mathbf{u}, \mathbf{u}}^{-1} \mathbf{u}, \text{diag}(\mathbf{K}_{\mathbf{f}, \mathbf{f}} - \mathbf{Q}_{\mathbf{f}, \mathbf{f}})) qFITC(f∗∣u)=p(f∗∣u)=N(K∗,uKu,u−1u,K∗,∗−Q∗,∗) \begin{aligned} q_{\text{FITC}}(\mathbf{f}_* | \mathbf{u}) &= p(\mathbf{f}_* | \mathbf{u}) \\ &= \mathcal{N}(\mathbf{K}_{*, \mathbf{u}} \mathbf{K}_{\mathbf{u}, \mathbf{u}}^{-1} \mathbf{u}, \mathbf{K}_{*, *} - \mathbf{Q}_{*, *}) \end{aligned}

となる。今までと同じように計算して、同時分布

qFITC(f,f∗)=N(0,(Qf,f+diag(Kf,f−Qf,f)Q∗,fQf,∗K∗,∗)) \begin{aligned} q_{\text{FITC}}(\mathbf{f}, \mathbf{f}_*) = \mathcal{N}\left( \mathbf{0}, \left( \begin{matrix} \mathbf{Q}_{\mathbf{f}, \mathbf{f}} + \text{diag}(\mathbf{K}_{\mathbf{f}, \mathbf{f}} - \mathbf{Q}_{\mathbf{f}, \mathbf{f}}) & \mathbf{Q}_{*, \mathbf{f}} \\ \mathbf{Q}_{\mathbf{f}, *} & \mathbf{K}_{*, *} \end{matrix} \right) \right) \end{aligned}

と予測分布

qFITC(f∗∣y)=N(Q∗,f(Qf,f+Λ)−1y,K∗,∗−Q∗,f(Qf,f+Λ)−1Qf,∗)=N(K∗,uΣKu,fΛ−1y,K∗,∗−Q∗,∗+K∗,uΣKu,∗) \begin{aligned} q_{\text{FITC}}(\mathbf{f}_* | \mathbf{y}) & = \mathcal{N}(\mathbf{Q}_{*, \mathbf{f}}(\mathbf{Q}_{\mathbf{f},\mathbf{f}} + \Lambda)^{-1} \mathbf{y}, \mathbf{K}_{*, *} - \mathbf{Q}_{*, \mathbf{f}}(\mathbf{Q}_{\mathbf{f},\mathbf{f}} + \Lambda)^{-1} \mathbf{Q}_{\mathbf{f}, *}) \\ & = \mathcal{N} (\mathbf{K}_{*, \mathbf{u}} \Sigma \mathbf{K}_{\mathbf{u}, \mathbf{f}} \Lambda^{-1} \mathbf{y}, \mathbf{K}_{*, *} - \mathbf{Q}_{*, *} + \mathbf{K}_{*, \mathbf{u}} \Sigma \mathbf{K}_{\mathbf{u}, *}) \end{aligned}

を得る。ここで、

Σ=(Ku,fΛ−1Kf,u+Ku,u)−1,Λ=diag(Kf,f−Qf,f+σ2I) \begin{aligned} \Sigma & = (\mathbf{K}_{\mathbf{u}, \mathbf{f}} \Lambda^{-1} \mathbf{K}_{\mathbf{f}, \mathbf{u}} + \mathbf{K}_{\mathbf{u}, \mathbf{u}})^{-1}, \\ \Lambda &= \text{diag}(\mathbf{K}_{\mathbf{f}, \mathbf{f}} - \mathbf{Q}_{\mathbf{f}, \mathbf{f}} + \sigma^2 \mathbf{I}) \end{aligned}

である。

最後に

統一的な枠組みで SoD, SoR, DTC, FITC と順に読んでいくことで補助変数法を一望することができて面白かった。

本文では、一部の相関を無視しないで考える The Partially Independent Training Conditional (PITC) や Inducing Variables の選び方についても触れられているようだ。

時間があれば、各補助変数法を実際に実装してみて、回帰結果がどの程度変わるのか確認してみたい。

-> やりました

ガウス過程の補助変数法をJuliaで実装、回帰結果を比較