我已经使用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) = 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
的值现在不下溢,但在某些情况下它仍然可以下溢。 det(c * A) = (c ^ n) * det(A)
,其中A
是n
矩阵的n
,然后将您的矩阵与c
预先相乘,计算行列式,然后除以c ^ n
获取实际的。 更新资料
我考虑了另一种选择。如果在#1或#2的最后阶段仍然经常发生下溢,那么最好专门针对这些最后的操作提高精度,例如,通过使用 GNU MPFR 。