行列式をLU分解で求める
EMアルゴリズムで多変量な混合正規分布のパラメータを推定するプログラムを書いてみようと思ったのだが、多変量正規分布の中に行列式を計算する箇所があり (以下の式の|Σ|)、 行列式の計算をどのように実装するのかわからなかったので調べて 実装してみた (gist)。
\[ \frac{1}{\sqrt{2\pi^D |\Sigma|}} \exp\{-\frac{1}{2}({\mathbf x}-{\mathbf \mu})^T\Sigma^{-1}({\mathbf x}-{\mathbf \mu})\} \]
NumPyのdet関数を覗いてみると"The determinant is computed via LU factorization using the LAPACK routine z/dgetrf.“とのこと。 LU分解を使って計算している。恥ずかしながら、LU分解は名前しか聞いたことが無かったが、行列を上三角行列Uと下三角行列Lに分解することらしい。 なぜLU分解で行列式を求めるのかを調べると、行列の積の行列式は、行列式の積に等しくて (see 行列の積の行列式)、三角行列の行列式は対角成分の積に等しいため。 計算の途中で、Uの対角成分が小さい場合にゼロ除算が起こりうるので、事前に行列の行を、対角成分が大きくなるように並び替える処理が必要。そうして求めたLとUは、並び替え後の行列に対するものなので、並び替える前のもとの行列の行列式を求める際には、並び替え後の行列の行列式に対して、並び替えた回数だけ-1をかける。
