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;

现在，在某些矩阵上，我不幸观察到数值下溢，因此行列式输出为零（实际上不是）。

我怎样才能避免这种下溢？

一种解决方案是获取行列式的对数，而不是行列式。然而，

  • 我不知道该怎么做；
  • 我怎么能把这些日志加起来？

任何帮助是极大的赞赏。

4

2 回答 2

4

我想到了两个主要选项：

  1. 方阵特征值的乘积就是这个矩阵的行列式，因此每个特征值的对数之和就是这个矩阵行列式的对数。假设det(A) = a和det(B) = b为紧凑符号。A在对 2 个矩阵和应用上述之后B，我们最终得到log(a)and 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), where Ais nby nmatrix ，并将矩阵与 some 预乘c，计算行列式，然后将其除以c ^ n得到实际的行列式.

更新

我又想了一个选择。如果在 #1 或 #2 的最后阶段您仍然过于频繁地遇到下溢，那么专门为这些最后的操作提高精度可能是个好主意，例如，通过使用GNU MPFR。

于 2013-12-11T21:16:39.897 回答
1

您可以使用 Householder 消除来获得 的 QR 分解delta0。那么 Q 部分的行列式是 +/-1（取决于您是否进行了偶数或奇数次反射），而 R 部分的行列式是对角线元素的乘积。这两个都很容易计算而不会陷入下溢地狱——你甚至可能不关心第一个。

于 2013-12-12T03:58:11.310 回答