c++ - 在 Eigen 中获取大矩阵的行列式时避免数值下溢

Question

我已经使用 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;

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

如何避免这种下溢？

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

• 我不知道该怎么做；
• 我怎样才能将这些日志相加？

非常感谢任何帮助。

Answer 1

我想到了两个主要选项:

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 或#2 的最后阶段您仍然经常遇到下溢，那么专门为这些最后的操作提高精度可能是个好主意，例如，通过利用 GNU MPFR .