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

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

如何避免这种下溢?

一种解决方案是获取行列式的对数,而不是行列式。然而,
  • 我不知道该怎么做;
  • 然后如何添加这些日志?

  • 任何帮助是极大的赞赏。

    最佳答案

    我想到了2个主要选项:

  • 方阵特征值的乘积是该矩阵的行列式,因此,每个特征值的对数之和是该矩阵行列式的对数。假定det(A) = adet(B) = b为紧凑表示法。在对2个矩阵AB应用上述条件后,我们得到log(a)log(b),那么实际上是正确的:
    log(a + b) = log(a) + log(1 + e ^ (log(b) - log(a)))
    

    是的,我们得到总和的对数。接下来您将如何处理?我不知道,取决于您必须执行的操作。如果必须通过e ^ log(a + b) = a + b删除对数,则可能很幸运a + b的值现在不下溢,但在某些情况下它仍然可以下溢。
  • 执行聪明的预处理;这里可能有很多选择,您最好从一些可信任的来源那里阅读它们,因为这是一个严肃的话题。针对该特定问题进行预处理的最简单(可能是最便宜的)示例可能是调用det(c * A) = (c ^ n) * det(A),其中An矩阵的n,然后将您的矩阵与c预先相乘,计算行列式,然后除以c ^ n获取实际的。

  • 更新资料

    我考虑了另一种选择。如果在#1或#2的最后阶段仍然经常发生下溢,那么最好专门针对这些最后的操作提高精度,例如,通过使用 GNU MPFR

    07-24 09:45
    查看更多