手帳と試行

学んだことをアウトプットしていきます。 日々、ノートあるのみ。
2023-01-30
回帰モデル

事後予測分布とWoodburyの公式

事後予測分布の式を、もう少し簡単にする。

目次
Woodburyの公式 事後予測分布の式の書き換え 1. 準備 2. 共分散行列の計算 3. 期待値の計算 結果 ちなみに

Woodburyの公式

逆行列の計算に便利な次の公式がある。

(A+BDC)−1=A−1−A−1B(D−1+CA−1B)CA−1\begin{aligned} (\bm A + \bm B \bm D \bm C)^{-1} = \bm A^{-1} - \bm A^{-1} \bm B (\bm D^{-1} + \bm C \bm A^{-1} \bm B) \bm C \bm A^{-1} \end{aligned}
A∈RN×NB∈RN×MC∈RM×ND∈RM×M\begin{aligned} \bm A \in{}& \R^{N \times N} \\ \bm B \in{}& \R^{N \times M} \\ \bm C \in{}& \R^{M \times N} \\ \bm D \in{}& \R^{M \times M} \end{aligned}

特に D=IM\bm D = \bm I_M の場合には次式に帰着する。

(A+BC)−1=A−1−A−1B(IM+CA−1B)CA−1\begin{aligned} (\bm A + \bm B \bm C)^{-1} = \bm A^{-1} - \bm A^{-1} \bm B (\bm I_M + \bm C \bm A^{-1} \bm B) \bm C \bm A^{-1} \end{aligned}

導出は省略する。

特に覚える必要はない。というか、使うときは徹底的に使いまくるので勝手に覚えることになる。

事後予測分布の式の書き換え

事後予測分布は次のようであった。

p(y∗∣X∗,X,y)=Nd∗(y∗∣my∗∣y,Vy∗∣y)\begin{aligned} p(\bm y_\ast | \bm X_\ast, \bm X, \bm y) = \mathcal N_{d^\ast} (\bm y_\ast | \bm m_{y_\ast | y}, \bm V_{y_\ast | y}) \end{aligned}
{my∗∣y=1σ2Vy∗∣yX∗V+Vd−1mdVy∗∣y−1=1σ2(Id∗−1σ2X∗V+X∗T)V+−1=1σ2X∗TX∗+Vd−1md=Vd(1σ2XTy+V0−1m0)Vd−1=1σ2XTX+V0−1\left\{\begin{aligned} \bm m_{y_\ast | y} ={}& \frac{1}{\sigma^2} \bm V_{y_\ast | y} \bm X_\ast \bm V_{+} \bm V_d^{-1} \bm m_d \\ \bm V_{y_\ast | y}^{-1} ={}& \frac{1}{\sigma^2} \left( \bm I_{d^\ast} - \frac{1}{\sigma^2} \bm X_\ast \bm V_+ \bm X_\ast^\mathsf{T} \right) \\ \bm V_+^{-1} ={}& \frac{1}{\sigma^2} \bm X_\ast^\mathsf{T} \bm X_\ast + \bm V_d^{-1} \\ \bm m_d ={}& \bm V_d \left( \frac{1}{\sigma^2} \bm X^\mathsf{T} \bm y + \bm V_0^{-1} \bm m_0 \right) \\ \bm V_d^{-1} ={}& \frac{1}{\sigma^2} \bm X^\mathsf{T} \bm X + \bm V_0^{-1} \end{aligned}\right.

出てくる記号が多くて計算が大変である。というか、逆行列計算が多すぎて、計算機を使うにしてもちょっと現実的ではない。

そこで、Woodburyの公式を適用してもう少し簡単な形に書き直していく。

1. 準備

まず Vd\bm V_d、md\bm m_d、さらに V+\bm V_+ を計算する。

Vd=(1σ2XTX+V0−1)−1=V0−V0XT(XV0XT+σ2Id)−1XV0md=Vd(1σ2XTy+V0−1m0)=(V0−V0XT(XV0XT+σ2Id)−1XV0)(1σ2XTy+V0−1m0)=1σ2V0XT(Id−(XV0XT+σ2Id)−1XV0XT)⏟(1)y+(In−V0XT(XV0XT+σ2Id)−1X)m0∣(1)=Id⏟A−1−(XV0XT⏟C+σ2Id⏟D−1)−1XV0XT⏟C=(Id⏟A+1σ2XV0XT⏟BDC)−1=σ2(XV0XT+σ2Id)−1=V0XT(XV0XT+σ2Id)−1y+(In−V0XT(XV0XT+σ2Id)−1X)m0V+=(1σ2X∗TX∗+Vd−1)−1=Vd−VdX∗T(X∗VdX∗T+σ2Id∗)−1X∗Vd\begin{aligned} \bm V_d ={}& \left( \frac{1}{\sigma^2} \bm X^\mathsf{T} \bm X + \bm V_0^{-1} \right)^{-1} \\ ={}& \bm V_0 - \bm V_0 \bm X^\mathsf{T} (\bm X \bm V_0 \bm X^\mathsf{T} + \sigma^2 \bm I_d)^{-1} \bm X \bm V_0 \\ \\ \bm m_d ={}& \bm V_d \left( \frac{1}{\sigma^2} \bm X^\mathsf{T} \bm y + \bm V_0^{-1} \bm m_0 \right) \\ ={}& \left( \bm V_0 - \bm V_0 \bm X^\mathsf{T} (\bm X \bm V_0 \bm X^\mathsf{T} + \sigma^2 \bm I_d)^{-1} \bm X \bm V_0 \right) \left( \frac{1}{\sigma^2} \bm X^\mathsf{T} \bm y + \bm V_0^{-1} \bm m_0 \right) \\ ={}& \frac{1}{\sigma^2} \bm V_0 \bm X^\mathsf{T} \underbrace{ \left( \bm I_d - (\bm X \bm V_0 \bm X^\mathsf{T} + \sigma^2 \bm I_d)^{-1} \bm X \bm V_0 \bm X^\mathsf{T} \right) }_{(1)} \bm y \\ & + \left( \bm I_n - \bm V_0 \bm X^\mathsf{T} (\bm X \bm V_0 \bm X^\mathsf{T} + \sigma^2 \bm I_d)^{-1} \bm X \right) \bm m_0 \\ &\left|\small\quad\begin{aligned} (1) ={}& \underbrace{\bm I_d}_{\bm A^{-1}} - ( \underbrace{ \bm X \bm V_0 \bm X^\mathsf{T} }_{ \bm C } + \underbrace{ \sigma^2 \bm I_d }_{\bm D^{-1}} )^{-1} \underbrace{ \bm X \bm V_0 \bm X^\mathsf{T} }_{ \bm C } \\ ={}& ( \underbrace{ \bm I_d }_{\bm A} + \underbrace{ \frac{1}{\sigma^2} \bm X \bm V_0 \bm X^\mathsf{T} }_{ \bm B \bm D \bm C} )^{-1} \\ ={}& \sigma^2 \left( \bm X \bm V_0 \bm X^\mathsf{T} + \sigma^2 \bm I_d \right)^{-1} \\ \end{aligned}\right. \\ ={}& \bm V_0 \bm X^\mathsf{T} \left( \bm X \bm V_0 \bm X^\mathsf{T} + \sigma^2 \bm I_d \right)^{-1} \bm y + \left( \bm I_n - \bm V_0 \bm X^\mathsf{T} (\bm X \bm V_0 \bm X^\mathsf{T} + \sigma^2 \bm I_d)^{-1} \bm X \right) \bm m_0 \\ \\ \bm V_+ ={}& \left( \frac{1}{\sigma^2} \bm X_\ast^\mathsf{T} \bm X_\ast + \bm V_d^{-1} \right)^{-1} \\ ={}& \bm V_d - \bm V_d \bm X_\ast^\mathsf{T} (\bm X_\ast \bm V_d \bm X_\ast^\mathsf{T} + \sigma^2 \bm I_{d^\ast})^{-1} \bm X_\ast \bm V_d \\ \end{aligned}
2. 共分散行列の計算

続いて Vy∗∣y\bm V_{y_\ast | y} を計算する。

Vy∗∣y=σ2(Id∗−1σ2X∗V+X∗T)−1=σ2(Id∗−X∗(−σ2V+−1+X∗TX∗)−1X∗T)=X∗(V+−1−1σ2X∗TX∗)−1⏟VdX∗T+σ2Id∗=X∗VdX∗T+σ2Id∗=X∗V0X∗T+σ2Id∗−X∗V0XT(XV0XT+σ2Id)−1XV0X∗T\begin{aligned} \bm V_{y_\ast | y} ={}& \sigma^2 \left( \bm I_{d^\ast} - \frac{1}{\sigma^2} \bm X_\ast \bm V_+\bm X_\ast^\mathsf{T} \right)^{-1} \\ ={}& \sigma^2 \left( \bm I_{d^\ast} - \bm X_\ast (-\sigma^2 \bm V_+^{-1} + \bm X_\ast^\mathsf{T} \bm X_\ast)^{-1} \bm X_\ast^\mathsf{T} \right) \\ ={}& \bm X_\ast \underbrace{\left( \bm V_+^{-1} - \frac{1}{\sigma^2} \bm X_\ast^\mathsf{T} \bm X_\ast \right)^{-1}}_{\bm V_d} \bm X_\ast^\mathsf{T} + \sigma^2 \bm I_{d^\ast} \\ ={}& \bm X_\ast \bm V_d \bm X_\ast^\mathsf{T} + \sigma^2 \bm I_{d^\ast} \\ ={}& \bm X_\ast \bm V_0 \bm X_\ast^\mathsf{T} + \sigma^2 \bm I_{d^\ast} - \bm X_\ast \bm V_0 \bm X^\mathsf{T} ( \bm X \bm V_0 \bm X^\mathsf{T} + \sigma^2 \bm I_d )^{-1} \bm X \bm V_0 \bm X_\ast^\mathsf{T} \end{aligned}
3. 期待値の計算

最後に my∗∣y\bm m_{y_\ast | y}。これはちょっと面倒くさい。

my∗∣y=1σ2Vy∗∣yX∗V+Vd−1md∣V+=Vd−VdX∗T(σ2Id∗+X∗VdX∗T⏟Vy∗∣y)−1X∗Vd=Vd−VdX∗TVy∗∣y−1X∗Vd=1σ2Vy∗∣yX∗(In−VdX∗TVy∗∣y−1X∗)md=1σ2Vy∗∣y(Id∗−X∗VdX∗TVy∗∣y−1)X∗md=1σ2Vy∗∣y(Vy∗∣y−X∗VdX∗T)⏟σ2Id∗Vy∗∣y−1X∗md=X∗md=X∗V0XT(XV0XT+σ2Id)−1y+(X∗−X∗V0XT(XV0XT+σ2Id)−1X)m0\begin{aligned} \bm m_{y_\ast | y} ={}& \frac{1}{\sigma^2} \bm V_{y_\ast | y} \bm X_\ast \bm V_+ \bm V_d^{-1} \bm m_d \\ &\left|\small\quad\begin{aligned} \bm V_+ ={}& \bm V_d - \bm V_d \bm X_\ast^\mathsf{T} (\underbrace{\sigma^2 \bm I_{d^\ast} + \bm X_\ast \bm V_d \bm X_\ast^\mathsf{T}}_{\bm V_{y_\ast | y}})^{-1} \bm X_\ast \bm V_d \\ ={}& \bm V_d - \bm V_d \bm X_\ast^\mathsf{T} \bm V_{y_\ast | y}^{-1} \bm X_\ast \bm V_d \end{aligned}\right. \\ ={}& \frac{1}{\sigma^2} \bm V_{y_\ast | y} \bm X_\ast (\bm I_n - \bm V_d \bm X_\ast^\mathsf{T} \bm V_{y_\ast | y}^{-1} \bm X_\ast) \bm m_d \\ ={}& \frac{1}{\sigma^2} \bm V_{y_\ast | y} (\bm I_{d^\ast} - \bm X_\ast \bm V_d \bm X_\ast^\mathsf{T} \bm V_{y_\ast | y}^{-1}) \bm X_\ast \bm m_d \\ ={}& \frac{1}{\sigma^2} \bm V_{y_\ast | y} \underbrace{(\bm V_{y_\ast | y} - \bm X_\ast \bm V_d \bm X_\ast^\mathsf{T})}_{\sigma^2 \bm I_{d^\ast}} \bm V_{y_\ast | y}^{-1}\bm X_\ast \bm m_d \\ ={}& \bm X_\ast \bm m_d \\ ={}& \bm X_\ast \bm V_0 \bm X^\mathsf{T} \left( \bm X \bm V_0 \bm X^\mathsf{T} + \sigma^2 \bm I_d \right)^{-1} \bm y \\ & + \left( \bm X_\ast - \bm X_\ast \bm V_0 \bm X^\mathsf{T} (\bm X \bm V_0 \bm X^\mathsf{T} + \sigma^2 \bm I_d)^{-1} \bm X \right) \bm m_0 \\ \end{aligned}

計算の途中で my∗∣y=X∗md\bm m_{y_\ast | y} = \bm X_\ast \bm m_d というものが出てきた。これは、y=Xw+ε\bm y = \bm X \bm w + \bm \varepsilon において、

y→my∗∣y,X→X∗,w→md\bm y \to \bm m_{y_\ast | y}, \quad\bm X \to \bm X_\ast, \quad\bm w \to \bm m_d

と書き換えたものに相当する。

このような置き換えはなかなか素朴なものに感じられるが、実際には「確率的モデル仮定→事後分布の導出→事後予測分布の導出」というステップをきちんと踏んだ上で得られた結果である。妥当性がそれなりに担保されているといえるだろう。

結果

ということで、事後予測分布は次のような形で書ける。

p(y∗∣X∗,X,y)=Nd∗(y∗∣my∗∣y,Vy∗∣y)\begin{aligned} p(\bm y_\ast | \bm X_\ast, \bm X, \bm y) = \mathcal N_{d^\ast} (\bm y_\ast | \bm m_{y_\ast | y}, \bm V_{y_\ast | y}) \end{aligned}
{my∗∣y=X∗V0XT(XV0XT+σ2Id)−1y+(X∗−X∗V0XT(XV0XT+σ2Id)−1X)m0Vy∗∣y=X∗V0X∗T+σ2Id∗−X∗V0XT(XV0XT+σ2Id)−1XV0X∗T\left\{\begin{aligned} \bm m_{y_\ast | y} ={}& \bm X_\ast \bm V_0 \bm X^\mathsf{T} \left( \bm X \bm V_0 \bm X^\mathsf{T} + \sigma^2 \bm I_d \right)^{-1} \bm y \\ & + \left( \bm X_\ast - \bm X_\ast \bm V_0 \bm X^\mathsf{T} (\bm X \bm V_0 \bm X^\mathsf{T} + \sigma^2 \bm I_d)^{-1} \bm X \right) \bm m_0 \\ \bm V_{y_\ast | y} ={}& \bm X_\ast \bm V_0 \bm X_\ast^\mathsf{T} + \sigma^2 \bm I_{d^\ast} - \bm X_\ast \bm V_0 \bm X^\mathsf{T} ( \bm X \bm V_0 \bm X^\mathsf{T} + \sigma^2 \bm I_d )^{-1} \bm X \bm V_0 \bm X_\ast^\mathsf{T} \\ \end{aligned}\right.

ちなみに

パラメータの事前分布の期待値が m0=0\bm m_0 = \bm 0 である場合、事後予測分布の期待値と共分散行列は

my∗∣y=X∗V0XT(XV0XT+σ2Id)−1yVy∗∣y=X∗V0X∗T+σ2Id∗−X∗V0XT(XV0XT+σ2Id)−1XV0X∗T\begin{aligned} \bm m_{y_\ast | y} ={}& \bm X_\ast \bm V_0 \bm X^\mathsf{T} (\bm X \bm V_0 \bm X^\mathsf{T} + \sigma^2 \bm I_d)^{-1} \bm y \\ \bm V_{y_\ast | y} ={}& \bm X_\ast \bm V_0 \bm X_\ast^\mathsf{T} + \sigma^2 \bm I_{d^\ast} - \bm X_\ast \bm V_0 \bm X^\mathsf{T} ( \bm X \bm V_0 \bm X^\mathsf{T} + \sigma^2 \bm I_d )^{-1} \bm X \bm V_0 \bm X_\ast^\mathsf{T} \\ \end{aligned}

となる。面白いことに、これはいくつかのかたまりに分けて解釈することができる。

my∗∣y=X∗V0XT⏟K∗T(XV0XT+σ2Id)−1⏟(K+σ2Id)−1yVy∗∣y=X∗V0X∗T+σ2Id∗⏟K∗∗+σ2Id∗−X∗V0XT⏟K∗T(XV0XT+σ2Id)−1⏟(K+σ2Id)−1XV0X∗T⏟K∗\begin{aligned} \bm m_{y_\ast | y} ={}& \underbrace{ \bm X_\ast \bm V_0 \bm X^\mathsf{T} }_{\bm K_\ast^\mathsf{T}} \underbrace{ (\bm X \bm V_0 \bm X^\mathsf{T} + \sigma^2 \bm I_d)^{-1} }_{(\bm K + \sigma^2 \bm I_d)^{-1}} \bm y \\ \bm V_{y_\ast | y} ={}& \underbrace{ \bm X_\ast \bm V_0 \bm X_\ast^\mathsf{T} + \sigma^2 \bm I_{d^\ast} }_{\bm K_{\ast\ast} + \bm \sigma^2 \bm I_{d^\ast}} - \underbrace{ \bm X_\ast \bm V_0 \bm X^\mathsf{T} }_{\bm K_\ast^\mathsf{T}} \underbrace{ ( \bm X \bm V_0 \bm X^\mathsf{T} + \sigma^2 \bm I_d )^{-1} }_{(\bm K + \sigma^2 \bm I_d)^{-1}} \underbrace{ \bm X \bm V_0 \bm X_\ast^\mathsf{T} }_{\bm K_\ast} \\ \end{aligned}

これらをあらかじめ計算しておけば、事後予測分布を

p(y∗∣X∗,X,y)=Nd∗(y∗∣my,Vy)p(\bm y_\ast | \bm X_\ast, \bm X, \bm y) = \mathcal N_{d^\ast} (\bm y_\ast | \bm m_y, \bm V_y)
my∗∣y=K∗T(K+σ2Id)−1yVy∗∣y=K∗∗+σ2Id∗−K∗T(K+σ2Id)−1K∗\begin{aligned} \bm m_{y_\ast | y} ={}& \bm K_\ast^\mathsf{T} (\bm K + \sigma^2 \bm I_d)^{-1} \bm y \\ \bm V_{y_\ast | y} ={}& \bm K_{\ast\ast} + \sigma^2 \bm I_{d^\ast} - \bm K_\ast^\mathsf{T} (\bm K + \sigma^2 \bm I_d)^{-1} \bm K_\ast \\ \end{aligned}
{K=XV0XTK∗=XV0X∗TK∗∗=X∗V0X∗T\left\{\begin{aligned} \bm K ={}& \bm X \bm V_0 \bm X^\mathsf{T} \\ \bm K_{\ast} ={}& \bm X \bm V_0 \bm X_\ast^\mathsf{T} \\ \bm K_{\ast\ast} ={}& \bm X_\ast \bm V_0 \bm X_\ast^\mathsf{T} \\ \end{aligned}\right.

と、比較的シンプルな形で記述することができる。Woodburyの公式を適用する前のものに比べてば、圧倒的に計算が簡単になっていることがわかるだろう。

しかし、それでも計算量が大きいという問題が立ちはだかる。ボトルネックは逆行列 K−1\bm K^{-1} の計算部分であり、時間計算量が O(d3)\mathcal O(d^3)、空間計算量が O(d2)\mathcal O(d^2) だけ要求される。そのためデータの個数 dd が大きくなると、現実的な時間では計算が終わらなくなるおそれがある。

回帰モデル