共分散行列 (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();

参考

既に実行中のジョブを走らせたまま端末を閉じる際の手順メモ

計算サーバ等で既に実行し始めてしまった時間のかかるプロセスについて、走らせ続けたまま端末を閉じたいことがある。

実行し始める前のプロセスについては、 nohup を使った方法がある。

ここでは、通常のコマンドで実行し始めてしまったプロセスについて考える。

disown を使う

時間がかかるタスクを通常のコマンドで実行し始めてしまったとする:

$ ./heavy_task.sh

対象のジョブを一旦停止する:

^Z

対象のジョブをバックグラウンド実行する:

$ bg

ジョブ番号を確認する:

$ jobs

disown コマンドによってバックグラウンド実行中のジョブをジョブテーブルから除外する:

$ disown %1

(これは対象のジョブ番号が1だったときの入力例)

以上の操作によって、端末を閉じてもプロセスは残ったままになる。

参考:

対数正規分布 (log-normal distribution) の可視化

0よりも大きい値しかとらない確率変数*1について、その事前分布を設定するときに、対数正規分布 (log-normal distribution) を用いることがある。

定義

対数正規分布とは、ある確率変数  x について、その対数  \log x の分布が正規分布に従うような分布のことである。以下の確率密度関数を持つ。

\displaystyle
\mathcal{LN}(x \mid \mu, \sigma^{2} ) = \frac{1}{x} \cdot \frac{1}{\sqrt{2 \pi \sigma^{2}}} \cdot \exp \left\{ - \frac{ ( \log x - \mu )^{2} }{2 \sigma^{2} } \right\}

ただし、  x \in \mathbb{R}_{> 0}, \mu \in \mathbb{R}, \sigma \in \mathbb{R}_{> 0} である。

性質

以下の性質を持つ。

  • meanは  \exp(\mu + \sigma^{2} / 2) で表される。
  • medianは  \exp(\mu) で表される
  • modeは  \exp(\mu - \sigma^{2}) で表される

直感的には、  \sigma^{2} が比較的小さいとき、mean、median、modeはいずれも  \exp(\mu) に近くなり、これを中心とした左右対称に近い分布になっていると思われる。逆に  \sigma^{2} を大きくしていくと、mean(つまり重心)は  \exp(\mu) よりも大きい方にずれ、mode(つまり最大をとる値)は  \exp(\mu) よりも小さい方にずれていく。

可視化

対数正規分布の形は少し想像が難しい。以下に可視化例を示す。

f:id:yuki-koyama:20190728094015p:plain
対数正規分布 (log-normal distribution) の可視化(  \mu を固定)

f:id:yuki-koyama:20190728094242p:plain
対数正規分布 (log-normal distribution) の可視化(  \sigma^{2} を固定)

以下に可視化に用いたPythonコードを示す。

おまけ:対数正規分布の対数

MAP推定などの統計的推定では、計算上対象とする分布の対数を扱うことが多い。対数正規分布の対数をとったものは

\displaystyle
\log \left\{ \mathcal{LN}(x \mid \mu, \sigma^{2} ) \right\} = - \log x - \frac{1}{2} \log(2 \pi \sigma^{2}) - \frac{ (\log x - \mu)^{2} }{2 \sigma^{2}}

となり、その微分

\displaystyle
\frac{d}{dx} \log \left\{ \mathcal{LN}(x \mid \mu, \sigma^{2} ) \right\} = - \frac{1}{x} \left(1 + \frac{\log x - \mu}{\sigma^{2}} \right)

となる。

*1:例えば分散  \sigma^{2} を確率変数とみなして統計的推定を行う場合など。

CMakeで管理するプロジェクトでEigen 3.Xを使うときのTips

Eigen 3.3以降

Eigen 3.3以降は

find_package(Eigen3 REQUIRED)

として、

target_link_libraries(mylibrary Eigen3::Eigen)

とすれば良い。

Eigen 3.3以降に関する公式ドキュメント: eigen.tuxfamily.org

なおmacOSのHomebrewで

brew install eigen

とすると、2019年3月15日現在、Eigen 3.3.7が入るので、この方法で良い。

Eigen 3.2以前

Eigen 3.2以前やEigen 3.3-betaなどは

find_package(Eigen3 REQUIRED)

としても、EIGEN3_INCLUDE_DIRなどの変数はセットされるが、上述のようなEigen3::Eigenというtargetは生成されない。そこで、

target_include_directories(mylibrary PUBLIC ${EIGEN3_INCLUDE_DIR})

などとしてヘッダーへのパスを通す必要がある。

なお、この方法だとmylibraryとEigenとの依存関係が直接的には指定できていないので、その後階層的に依存関係を構築したりする際にもしかするとどこかで弊害があるかもしれない(要確認)し、見た目にも依存関係が分かりにくい(主観)。

なおUbuntuのaptで

sudo apt install libeigen3-dev

とすると、Eigen 3.2やEigen 3.3-betaがインストールされる環境が存在するので、その場合はこのような方法をとる必要がある。

両方のケースに対応するCMake

両者に対応し、依存関係も直接指定するには、

find_package(Eigen3 REQUIRED)
if((NOT TARGET Eigen3::Eigen) AND (DEFINED EIGEN3_INCLUDE_DIR))
    add_library(AliasEigen3 INTERFACE)
    target_include_directories(AliasEigen3 INTERFACE ${EIGEN3_INCLUDE_DIR})
    add_library(Eigen3::Eigen ALIAS AliasEigen3)
endif()

とした上で

target_link_libraries(mylibrary Eigen3::Eigen)

とすれば良い。

Eigenはheader-only libraryなので、INTERFACEキーワードでtargetを設定した上で、ALIASキーワードを使ってEigen3::Eigenで指せるようにしている。

cmake.org

参考レポジトリ

下記のレポジトリの2019年3月15日現在の最新版では、上のような方法で両者のケースに対応している。

github.com github.com