概要
本記事では、コレスキー分解を用いて適切な変数変換を行うことで、多次元正規分布の確率密度関数に含まれる分散共分散行列の逆行列の数値計算を回避する方法について述べる。
多次元正規分布
\(M\) 次元正規分布の確率密度関数
$$\mathcal{N}(\mathbf{x}|\boldsymbol{\mu}, \boldsymbol{\Sigma})=\frac{1}{(2\pi)^\frac{M}{2}|\boldsymbol{\Sigma}|^\frac{1}{2}}\exp\left\{-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^\mathsf{T}\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})\right\}$$
の定義式には、分散共分散行列の逆行列 \(\boldsymbol{\Sigma}^{-1}\) の計算が含まれている。数値計算を行う際、速度や安定性の観点から、これを素直に計算するのは好ましくない。
指数関数の内部
平均・分散共分散行列・平均
$$(\mathbf{x}-\boldsymbol{\mu})^\mathsf{T}\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})\tag{1}$$
逆行列の計算を含む指数関数の内部は、多次元正規分布から派生する種々の計算でよく見る形状である。これを効率的に数値計算する方法を考える。
変数の置き換え
$$\mathbf{x}-\boldsymbol{\mu}=\mathbf{Az}\tag{2}$$
という変数の置き換えを行うと
$$(\mathbf{x}-\boldsymbol{\mu})^\mathsf{T}\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})=\mathbf{z}^\mathsf{T}\mathbf{A}^\mathsf{T}\boldsymbol{\Sigma}^{-1}\mathbf{Az}$$
と書ける。ここで、
$$\mathbf{A}^\mathsf{T}\boldsymbol{\Sigma}^{-1}\mathbf{A}=\mathbf{I}$$
という関係をみたす行列 \(A\) が存在することを仮定する（ここで、 \(\mathbf{I}\) は \(M\times M\) の単位行列である）。これを変形すると
$$\boldsymbol{\Sigma}^{-1}=(\mathbf{A}^\mathsf{T})^{-1}\mathbf{A}^{-1}$$
$$\boldsymbol{\Sigma}=\left\{(\mathbf{A}^\mathsf{T})^{-1}\mathbf{A}^{-1}\right\}^{-1}=\mathbf{A}\mathbf{A}^\mathsf{T}\tag{3}$$
という関係を意味することがわかる。
コレスキー分解
（特定の次元について分散が0とならない限り、）分散共分散行列は正定値対称行列であるので、以下のコレスキー分解を行うことができる。
$$\boldsymbol{\Sigma}=\mathbf{A}\mathbf{A}^\mathsf{T}$$
ここで、 \(\mathbf{A}\) は \(M\times M\) の下三角行列である。
※正式には分散共分散行列は半正定値対称行列であり、加えて浮動小数点数の誤差により、（半）正定値の条件をみたさなくなる場合がある。この場合にコレスキー分解を行うことはできないが、それに際しての対策はまた後日。
以上より、分散共分散行列に対して \((3)\) の関係をみたす行列 \(\mathbf{A}\) は存在する。
ベクトル \(z\) の導出
\((2)\) の関係式は
$$\mathbf{z}=\mathbf{A}^{-1}(\mathbf{x}-\boldsymbol{\mu})\tag{4}$$
と変形できる。 \(\mathbf{A}\) はコレスキー分解の結果求まる下三角行列であるため、この逆行列の計算は、分散共分散行列の逆行列の計算よりも簡単である。
下三角行列の逆行列
下三角行列の逆行列も下三角行列であり、その対角要素は、もとの対角要素の逆数になる。
（参考）
これを踏まえて、逆行列の定義式
$$\mathbf{A}\mathbf{A}^{-1}=\mathbf{I}$$
の要素を書き下すと
$$\begin{bmatrix}
a_{1,1} & 0 & \cdots & 0 & 0 & 0 \\
a_{2,1} & a_{2,2} & 0 & \cdots & 0 & 0 \\
a_{3,1} & a_{3,2} & a_{3,3} & 0 & \cdots & 0 \\
\vdots & & & & & \vdots \\
a_{M,1} & a_{M,2} & a_{M,3} & a_{M,4} & \cdots & a_{M,M}
\end{bmatrix}
\begin{bmatrix}
\frac{1}{a_{1,1}} & 0 & \cdots & 0 & 0 & 0 \\
b_{2,1} & \frac{1}{a_{2,2}} & 0 & \cdots & 0 & 0 \\
b_{3,1} & b_{3,2} & \frac{1}{a_{3,3}} & 0 & \cdots & 0 \\
\vdots & & & & & \vdots \\
b_{M,1} & b_{M,2} & b_{M,3} & b_{M,4} & \cdots & \frac{1}{a_{M,M}}
\end{bmatrix}=
\begin{bmatrix}
1 & 0 & \cdots & 0 & 0 & 0 \\
0 & 1 & 0 & \cdots & 0 & 0 \\
0 & 0 & 1 & 0 & \cdots & 0 \\
\vdots & & & & & \vdots \\
0 & 0 & 0 & 0 & \cdots & 1
\end{bmatrix}$$
より、逆行列の要素 \(b_{\cdot,\cdot}\) を、以下のように順次求めることができる。
$$\frac{a_{2,1}}{a_{1,1}}+a_{2,2}b_{2,1}=0\Leftrightarrow b_{2,1}=-\frac{1}{a_{2,2}}\frac{a_{2,1}}{a_{1,1}}$$
$$\frac{a_{3,1}}{a_{1,1}}+a_{3,2}b_{2,1}+a_{3,3}b_{3,1}=0\Leftrightarrow b_{3,1}=-\frac{1}{a_{3,3}}\left(\frac{a_{3,1}}{a_{1,1}}+a_{3,2}b_{2,1}\right)$$
$$\frac{a_{3,2}}{a_{2,2}}+a_{3,3}b_{3,2}=0\Leftrightarrow b_{3,2}=-\frac{1}{a_{3,3}}\frac{a_{3,2}}{a_{2,2}}$$
$$\vdots$$
ベクトル \(z\) の代入
\((4)\) を \((1)\) に代入すると、
$$(\mathbf{x}-\boldsymbol{\mu})^\mathsf{T}\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})=\mathbf{z}^\mathsf{T}\mathbf{A}^\mathsf{T}\boldsymbol{\Sigma}^{-1}\mathbf{A}\mathbf{z}=\mathbf{z}^\mathsf{T}\mathbf{z}$$
となる。途中で \((3)\) の関係式を用いた。
したがって、分散共分散行列の逆行列の計算は
- コレスキー分解
- 下三角行列の逆行列の計算
に置き換えて回避することができる。
確率密度関数全体
\(M\) 次元確率密度関数を
$$f(\mathbf{x})=\frac{1}{(2\pi)^\frac{M}{2}|\boldsymbol{\Sigma}|^\frac{1}{2}}\exp\left\{-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^\mathsf{T}\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})\right\}$$
とおく。確率密度関数において、確率変数を変換する際の法則に注意すると、
（参考）
今回は \((4)\) より
$$\mathbf{z}=\mathbf{A}^{-1}(\mathbf{x}-\boldsymbol{\mu})\Leftrightarrow\mathbf{x}=\mathbf{Az}+\boldsymbol{\mu}$$
という変換を行うので、変換後の確率密度関数 \(g(\mathbf{z})\) は
$$g(\mathbf{z})=|\mathbf{A}|f(\mathbf{Az}+\boldsymbol{\mu})$$
と書ける。これを計算して
$$g(\mathbf{z})=\frac{|\mathbf{A}|}{(2\pi)^\frac{M}{2}|\boldsymbol{\Sigma}|^\frac{1}{2}}\exp\left(-\frac{1}{2}\mathbf{z}^\mathsf{T}\mathbf{A}^\mathsf{T}\boldsymbol{\Sigma}^{-1}\mathbf{Az}\right)$$
$$=\frac{|\mathbf{A}|}{(2