2

Eigen ライブラリを使用して C++ で MCMC アルゴリズムを実装しました。アルゴリズムの主要部分は、最初にいくつかの行列計算が実行された後、結果の行列の行列式が取得され、出力に追加されるループです。例えば：

MatrixXd delta0;
NumericVector out(3);

out[0] = 0;
out[1] = 0;

for (int i = 0; i < s; i++) {
  ...
  delta0 = V*(A.cast<double>()-(A+B).cast<double>()*theta.asDiagonal());
  ...
  I = delta0.determinant()
  out[1] += I;
  out[2] += std::sqrt(I);
}
return out;

現在、特定の行列では、残念ながら数値のアンダーフローが観察されるため、行列式がゼロとして出力されます（実際にはそうではありません）。

このアンダーフローを回避するにはどうすればよいですか?

1 つの解決策は、行列式の代わりに、行列式のログを取得することです。でも、

  • これを行う方法がわかりません。
  • これらのログを追加するにはどうすればよいですか?

どんな助けでも大歓迎です。

4

2 に答える 2

4

私の頭に浮かぶ主なオプションは2つあります。

  1. 正方行列の固有値の積は、この行列の行列式であるため、各固有値の対数の和は、この行列の行列式の対数です。コンパクトな表記にはdet(A) = aとを仮定します。det(B) = b前述の 2 つの行列Aとを適用するとB、最終的にlog(a)とlog(b)となり、実際には次のことが成り立ちます。

    log(a + b) = log(a) + log(1 + e ^ (log(b) - log(a)))

    はい、合計の対数を取得します。次にそれをどうしますか？私は知りません、あなたがしなければならないことに依存します。で対数を削除する必要がある場合、 の値が現在アンダーフローしていないe ^ log(a + b) = a + bことは幸運かもしれませんが、場合によってはアンダーフローすることもあります。a + b

  2. 巧妙な前処理を実行します。ここにはたくさんのオプションがあるかもしれません。これは深刻なトピックであるため、信頼できる情報源からそれらについて読んだほうがよいでしょう。この特定の問題に対する前処理の最も単純な (そしておそらくこれまでで最も安価な) 例は、 をdet(c * A) = (c ^ n) * det(A)、ここAでnをn行列で思い出し、行列に を前もって掛け、行列式をc計算し、それを で割ってc ^ n実際の行列を取得することです。 .

アップデート

もう1つの選択肢を考えました。#1 または #2 の最後の段階でまだアンダーフローが頻繁に発生する場合は、たとえばGNU MPFRを使用して、これらの最後の操作の精度を特に高めることをお勧めします。

于 2013-12-11T21:16:39.897 に答える