共分散行列 (Covariance Matrix) の行列式 (determinant) の対数 (log) を計算する

f:id:yuki-koyama:20190815150220p:plain

下記ライブラリのガウス過程回帰 (Gaussian Process Regression) の実装の際に悩んだ点をメモしておく。

github.com

いつ行列式の対数 (log-determinant) を計算する必要があるか

例えば、多変量正規分布 (multivariate normal distribution) が関連する統計的推定において、対数尤度を計算する際などに、共分散行列の行列式の対数を計算する必要が生じる。

行列式 (determinant) を計算する際の問題点

特にサイズが大きな行列の行列式は大変小さな値をとることがあり、計算機の浮動小数点数の表現力がシビアに効いてくることも多い。これの対数をとると、計算が不安定になったり、infが問題となったりする。

正定値対称行列 (symmetric positive-definite matrix) の行列式の対数 (log-determinant) を直接計算する方法

共分散行列など、計算の対象とする行列が正定値対称行列である場合、行列式を明示的に計算せずに直接行列式の対数を計算することができる。正定値対称行列  \mathbf{K} をコレスキー分解によって

\displaystyle
\mathbf{K} = \mathbf{L}\mathbf{L}^{\textsf{T}}

と分解すると、  \mathbf{K} 行列式の対数は

\displaystyle
\log ( \det ( \mathbf{K} ) ) = 2 \cdot \text{sum} ( \log ( \text{diag2vec} ( \mathbf{L} ) ) )

によって計算することができる。ただし  \text{diag2vec} は対角要素を取り出してベクトルにする変換、右辺の  \log は要素毎の対数、  \text{sum} はベクトルの各要素の和を表す。(この式が成立することの証明は こちら を参照。)

C++でEigenを用いる場合の実装

ナイーブな方法(このようにするべきではない):

// #include <Eigen/LU>
// #include <cmath>

// const Eigen::MatrixXd K = ...
const double log_det_K = std::log(K.determinant());

工夫した方法(こちらの方が好ましいと思われる):

// #include <Eigen/Cholesky>

// const Eigen::MatrixXd K = ...
const double log_det_K = 2.0 * Eigen::LLT<Eigen::MatrixXd>(K).matrixL().toDenseMatrix().diagonal().array().log().sum();

参考