在特征值中获取大矩阵行列式时避免数值下溢

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。