多変量ガウス分布の周辺化は、かの有名なPRMLでも登場する概念です。
計算的にも結構しんどい内容が多く、意外と多くの人が飛ばしてしまっている内容なんじゃないのかなと思います。しかし、ガウス過程やカルマンフィルタの導出や概念を理解するには、この多変量ガウス分布の周辺化の理解は避けて通れません。
今回はこの多変量ガウス分布の周辺化を絵や図で説明・解説してみようと思います。
おそらく、この記事が多変量ガウス分布の周辺化について、世界一丁寧に解説している記事になると思います。
正直、この辺りは数学的にも難しく、概念的にもイメージしにくい内容となっています。この記事ではイメージしにくい、多変量ガウス分布の周辺化をできるだけ絵や図を用いて解説しており、ここまで詳しく解説している記事はないと思うので、PRML等の本で躓いた方はぜひ読んでみてください。
ちなみに、確率変数における周辺化について理解が怪しい方は、先にこちらの記事で確率変数の周辺化について解説しているので、先に読むことをお勧めします。
多変量ガウス分布の周辺化
まず、多変量ガウス分布の周辺化でやりたい内容を整理しましょう。
まず多変量ガウス分布の定義を下に示します。
\begin{equation}
\mathcal{N}(\bm{x}|\bm{μ}, \bm{\Sigma}) =
\frac{1}{\sqrt{(2\pi)^{D} |\Sigma|}}
exp
\biggl\{
-\frac{1}{2}(\bm{x}-\bm{μ})^T\Sigma^{-1}(\bm{x}-\bm{μ})
\biggr\}
\end{equation}ここらへんの内容は、下記の記事で詳しく解説しているので、ガウス分布について怪しい人はこちらについても参考にしてみてください。
今回は、$D$次元の多変量ガウス分布を考えることにします
数式の準備
今、多変量ガウス分布の定義である(1)の入力$\bm{x} \in \mathbb{R}$を次のように、$\bm{x}_1$と$\bm{x}_2$に分割することを考えます。
\begin{split}
\bm{x}_1 = (x_1, x_2, \dots, x_L) \\
\bm{x}_2= (x_{l+1}, \dots, x_D)
\end{split}イメージはこんな感じです。$\bm{x_1}$の次元数を$L$とすると、$\bm{x_2}$の次元数は$D-L$で表現できます。
また、(1)式における、平均ベクトル$\mu$と分散共分散行列$\Sigma$も、$\bm{x}$のように分割することを考えます。それぞれ次のようになります。
\begin{split}
\bm{\mu}_1 = (\mu_1, \mu_2, \dots, \mu_L) \\
\bm{\mu}_2= (\mu_{l+1}, \dots, \mu_D)
\end{split}続いて、分散共分散行列$\Sigma$も分解することを試みます。
$\bm{x}$や$\bm{\mu}$のように、
$\Sigma_{11}$、$\Sigma_{12}$、$\Sigma_{21}、$\Sigma_{22}$に分解します。
\begin{equation}
\Sigma =
\begin{pmatrix}
\sigma_{1^2} & \sigma_{12} & \dots & \sigma_{1L} & \sigma_{1 L+1} & \dots & \sigma_{1D}\\
\sigma_{21} & \sigma_{2^2} & \dots & \sigma_{2L} & \sigma_{2 L+1} & \dots & \sigma_{2D}\\
\vdots & \vdots & & \ddots & \vdots & & \vdots \\
\sigma_{L1} & \sigma_{L2} & \dots & \sigma_{L^2} & \sigma_{L L+1} & \dots & \sigma_{LD}\\
\sigma_{L+1~1} & \sigma_{L+1~2} & \dots & \sigma_{L+1~L} & \sigma_{{L+1}^2} & \dots & \sigma_{L+1 ~D}\\
\sigma_{L+2~1} & \sigma_{L+2~2} & \dots & \sigma_{L+2~L} & \sigma_{L+2~L+1} & \dots & \sigma_{L+2 ~D}\\
\vdots & \vdots & & \ddots & \vdots & & \vdots \\
\sigma_{D1} & \sigma_{D2} & \dots & \sigma_{DL} & \sigma_{D{L+1}} & \dots & \sigma_{D D}\\
\end{pmatrix}
\end{equation}を次のように分解します。
\begin{equation}
\begin{split}
\Sigma =
\begin{pmatrix}
\Sigma_{11} & \Sigma_{12} \\
\Sigma_{21} & \Sigma_{22} \\
\end{pmatrix}
\end{split}
\end{equation}$\Sigma_{11}, \Sigma_{12}$などがいきなり登場して、面食らってしまうかもしれませんが、単純に下記のように分割したと思ってください。
これは便宜的に、このように分割していることに注意してください。ちなみに、ガウス分布のパラメータである$\Sigma$は対称行列のため、次の等式が成り立つことに注意してください。
\begin{equation}
\begin{split}
\Sigma &= \Sigma^T \\
\Sigma_{12} &= \Sigma_{21}^T \\
\Sigma_{21} &= \Sigma_{12}^T
\end{split}
\end{equation}さらに、$\Sigma_{11}$と$\Sigma_{22}$も対称行列であることがわかります。(これのイメージがつきにくい人は、(2)の式をよく眺めてみてください。)
ここまででやっと数式の準備が整いました。
多変量ガウス分布の周辺化の公式
多変量ガウス分布の公式と仰々しく書きましたが、多変量ガウス分布の同時確率分布$p(\bm{x}) = p(\bm{x_1}, \bm{x_2})$ を、$p(\bm{x_2})$で周辺化すると、周辺分布$p(\bm{x_1})$は次のようになります。
\begin{equation}
\begin{split}
p(\bm{x_1}) &= \int p(\bm{x_1}, \bm{x_2}) d \bm{x_2} \\
&= \mathcal{N}(\bm{\mu_1, \Sigma_{11}})
\end{split}
\end{equation}多変量ガウス分布の周辺化の証明
では、多変量ガウス分布の周辺化の証明をしていきましょう。
まず、共分散行列$\Sigma$の逆行列を$\Lambda$と置きます。
\begin{equation}
\Lambda \equiv \Sigma^{-1}
\end{equation}今日分散行列の逆行列として定義した、$\Lambda$のことを精度行列と呼ぶことがあります。
共分散行列さえいまいち理解しにくいのに、さらに逆関数である精度行列にするなんて、なぜ？？と思う人がいるかもしれませんが、(1)式のがうす分布の定義式を見ると、精度行列の形でしか$\Sigma$が登場しないことがわかります。
とはいえ、少しわかりにくいので、次のように式変形するとわかりやすいでしょう。
\begin{equation}
\begin{split}
\mathcal{N}(\bm{x}|\bm{μ}, \bm{\Sigma}) &=
\frac{1}{\sqrt{(2\pi)^{D} |\Sigma|}}
exp
\biggl\{
-\frac{1}{2}(\bm{x}-\bm{μ})^T\Sigma^{-1}(\bm{x}-\bm{μ})
\biggr\} \\
&=
\frac{1}{\sqrt{(2\pi)^{D}}} |\Sigma|^{-\frac{1}{2}}
exp
\biggl\{
-\frac{1}{2}(\bm{x}-\bm{μ})^T\Sigma^{-1}(\bm{x}-\bm{μ})
\biggr\} \\
&=
\frac{1}{\sqrt{(2\pi)^{D}}} \Lambda^{-\frac{1}{2}}
exp
\biggl\{
-\frac{1}{2}(\bm{x}-\bm{μ})^T\Lambda(\bm{x}-\bm{μ})
\biggr\} \\
\end{split}
\end{equation}(7)式を見ると、$\Sigma$が$\Lambda$に置き換えられることがわかります。
このように、他変量ガウス分布の共分散行列を考える際には、よく精度行列を用いることがあるので、頭の片隅に入れておくと良いでしょう。
ここで、(6)式において、$\Sigma$が対称行列であることから、その逆行列である精度行列$\Lambda$も対称行列となります。
これは対称行列の有する基本的な性質なので、ここが怪しい人は下記の記事を参考にしてください。
ここで、精度行列$\Lambda$においても、(3)と同様な分解をすることを考えます。急にこのような分解が登場してしまいますが、以降の計算で必要のためなので、今は天下り的に理解すれば大丈夫です。
\begin{equation}
\begin{split}
\Lambda =
\begin{pmatrix}
\Lambda_{11} & \Lambda_{12} \\
\Lambda_{21} & \Lambda_{22} \\
\end{pmatrix}
\end{split}
\end{equation}このように分解します。
ここで$\Lambda$は対称行列なので、$\Lambda_{11}$、$\Lambda_{22}$は対称行列であり、$\Lambda_{12}^T = \Lambda_{21}$が成り立ちます。
ここで、注意してほしいのが、$\Sigma$と$\Lambda$の間には、$\Sigma$の逆行列が$\Lambda$であるという逆行列の関係性が成り立つのですが、分割した部分行列において、例えば$\Sigma_{11}$の逆行列が$\Lambda_{11}$になると言った関係は成り立たないことに注意してください。
ここらへんは、意識して次に進んでいきましょう。
ここで、精度行列に置き換えた多変量ガウス分布の指数部分に注目していきましょう。指数分布はこのように表現できています。
-\frac{1}{2}(\bm{x}-\bm{μ})^T\Lambda(\bm{x}-\bm{μ})ここで、
\bm{\mu} =
\begin{pmatrix}
\bm{\mu_1} \\
\bm{\mu_2} \\
\end{pmatrix},
\bm{x} =
\begin{pmatrix}
\bm{x_1} \\
\bm{x_2} \\
\end{pmatrix}であることから、上式と(8)を代入すると、
\begin{equation}
\begin{split}
-\frac{1}{2} &
\begin{pmatrix}
\bm{x_1} - \bm{\mu_1} \\
\bm{x_2} - \bm{\mu_2} \\
\end{pmatrix}^T
\begin{pmatrix}
\Lambda_{11} & \Lambda_{12} \\
\Lambda_{21} & \Lambda_{22} \\
\end{pmatrix}
\begin{pmatrix}
\bm{x_1} - \bm{\mu_1} \\
\bm{x_2} - \bm{\mu_2} \\
\end{pmatrix} \\
&=
-\frac{1}{2} ( \bm{x_1} - \bm{\mu_1})^T \Lambda_{11} ( \bm{x_1} - \bm{\mu_1})
-\frac{1}{2} ( \bm{x_1} - \bm{\mu_1})^T \Lambda_{12} ( \bm{x_2} - \bm{\mu_2}) \\
& ~~~~ -\frac{1}{2} ( \bm{x_2} - \bm{\mu_2})^T \Lambda_{21} ( \bm{x_1} - \bm{\mu_1})
-\frac{1}{2} ( \bm{x_2} - \bm{\mu_2})^T \Lambda_{22} ( \bm{x_2} - \bm{\mu_2}) \\
\end{split}
\end{equation}かなり煩雑な式であるが、多変量ガウス分布の確率密度関数(7)式の指数分は、分割した精度行列を用いてこのように表現きます。
さらに、ここを展開していきます。
この式展開では、 転置行列の和に関する線形性の性質を利用しています。
$A$, $B$を行列とするとき、
(A + B) ^T = A^T + B^T
詳しくは、こちらの記事をご覧ください。
また、対称行列に対して成り立つ、以下の演算も前提としています。
ベクトル$\bm{x}, \bm{y}$、対象行列$A$に対して、下記の交換則が成り立つ。
\begin{equation}
\bm{x}^T A \bm{y} = \bm{y}^T A\bm{x}
\end{equation}行列$A$は対称行列であることから、
A^T = A
が成り立つ。ここで、
\begin{split}
\bm{x}^T A \bm{y} &= \bm{x}^T (A \bm{y}) \\
&= (A \bm{y})^T \bm{x} \\
&= \bm{y}^T A^T \bm{x} \\
&= \bm{y}^T A \bm{x}
\end{split}以降の式変形では、これらの演算を前提として式変形を進めていきます。
(9)式をさらに展開すると、
\begin{equation}
\begin{split}
-\frac{1}{2} &( \bm{x_1} - \bm{\mu_1})^T \Lambda_{11} ( \bm{x_1} - \bm{\mu_1})
-\frac{1}{2} ( \bm{x_1} - \bm{\mu_1})^T \Lambda_{12} ( \bm{x_2} - \bm{\mu_2}) \\
& ~~~~ -\frac{1}{2} ( \bm{x_2} - \bm{\mu_2})^T \Lambda_{21} ( \bm{x_1} - \bm{\mu_1})
-\frac{1}{2} ( \bm{x_2} - \bm{\mu_2})^T \Lambda_{22} ( \bm{x_2} - \bm{\mu_2}) \\
&=
-\frac{1}{2} ( \bm{x_1}^T \Lambda_{11} \bm{x_1} - \bm{x_1}^T \Lambda_{11} \bm{\mu_1} -
\bm{\mu_1}^T \Lambda_{11} \bm{x_1} + \bm{\mu_1}^T \Lambda_{11} \bm{\mu_1} ) \\
& ~~~~ -\frac{1}{2} ( \bm{x_1}^T \Lambda_{12} \bm{x_2} - \bm{x_1}^T \Lambda_{12} \bm{\mu_2} - \bm{\mu_1}^T \Lambda_{12} \bm{x_2} + \bm{\mu_1}^T \Lambda_{12} \bm{\mu_2} ) \\
& ~~~~ -\frac{1}{2} ( \bm{x_2}^T \Lambda_{21} \bm{x_1} - \bm{x_2}^T \Lambda_{21} \bm{\mu_1} - \bm{\mu_2}^T \Lambda_{21} \bm{x_1} + \bm{\mu_2}^T \Lambda_{21} \bm{\mu_1} ) \\
& ~~~~ -\frac{1}{2} ( \bm{x_2}^T \Lambda_{22} \bm{x_2} - \bm{x_2}^T \Lambda_{22} \bm{\mu_2} - \bm{\mu_2}^T \Lambda_{22} \bm{x_2} + \bm{\mu_2}^T \Lambda_{22} \bm{\mu_2} ) \\
&= -\frac{1}{2} \bm{x_1}^T \Lambda_{11} \bm{x_1} + \bm{x_1}^T \Lambda_{11} \bm{\mu_1}
-\frac{1}{2} \bm{\mu_1}^T \Lambda_{11} \bm{\mu_1} \\
&~~~~ - \bm{x_2}^T \Lambda_{21} \bm{x_1} + \bm{x_2}^T \Lambda_{21} \bm{\mu_1} + \bm{\mu_2}^T \Lambda_{21} \bm{x_1} - \bm{\mu_2}^T \Lambda_{21} \bm{\mu_1} \\
&~~~~ -\frac{1}{2} \bm{x_2}^T \Lambda_{22} \bm{x_2} + \bm{x_2}^T \Lambda_{22} \bm{\mu_2}
-\frac{1}{2} \bm{\mu_2}^T \Lambda_{22} \bm{\mu_2}
\end{split}
\end{equation}ここまで、式変形をすることができました。
今回やりたいことは、$\bm{x}_2$における周辺分布を求めることです。周辺分布を求めるためには、$\bm{x}_2$について積分をする必要があるため、(11)式を$\bm{x}_2$について整理して、積分を行いやすいようにします。
まず、(11)式において、$\bm{x}_2$に関する項を抽出すると、
\begin{equation}
\begin{split}
& -\frac{1}{2} \bm{x_1}^T \Lambda_{11} \bm{x_1} + \bm{x_1}^T \Lambda_{11} \bm{\mu_1}
-\frac{1}{2} \bm{\mu_1}^T \Lambda_{11} \bm{\mu_1} \\
&~~~~ - \bm{x_2}^T \Lambda_{21} \bm{x_1} + \bm{x_2}^T \Lambda_{21} \bm{\mu_1} + \bm{\mu_2}^T \Lambda_{21} \bm{x_1} - \bm{\mu_2}^T \Lambda_{21} \bm{\mu_1} \\
&~~~~ -\frac{1}{2} \bm{x_2}^T \Lambda_{22} \bm{x_2} + \bm{x_2}^T \Lambda_{22} \bm{\mu_2}
-\frac{1}{2} \bm{\mu_2}^T \Lambda_{22} \bm{\mu_2} \\
&=
-\frac{1}{2} \bm{x_2}^T \Lambda_{22} \bm{x_2} + \bm{x_2}^T \Lambda_{22} \bm{\mu_2}
- \bm{x_2}^T \Lambda_{21} \bm{x_1} - \bm{x_2}^T \Lambda_{21} \bm{\mu_1} \\
&=
-\frac{1}{2} \bm{x_2}^T \Lambda_{22} \bm{x_2} + \bm{x_2}^T
\biggl \{
\Lambda_{22} \bm{\mu_2}
- \Lambda_{21} (\bm{x_1} - \bm{\mu_1})
\biggr \} \\
&=
-\frac{1}{2} \bm{x_2}^T \Lambda_{22} \bm{x_2} + \bm{x_2}^T
\bm{m}
\end{split}
\end{equation}と表現することができます。
ここで、(12)式の最後の式変形においては、$\bm{x}_2$の2次の項と1次の項にわけ、1次の項の係数を$\bm{m}$と置きました。
\begin{equation}
\bm{m} = \Lambda_{22} \bm{\mu_2}
- \Lambda_{21} (\bm{x_1} - \bm{\mu_1})
\end{equation}