CODEWITHSUNDEEP
c++floating-pointeigenunderflowdeterminants
Henrik
Henrik

Reputation: 14450

Avoid numerical underflow when obtaining determinant of large matrix in Eigen

I have implemented a MCMC algorithm in C++ using the Eigen library. The main part of the algorithm is a loop in which first some some matrix calculations are performed after which the determinant of the resulting matrix is obtained and added to the output. E.g.:

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;

Now on certain matrices I unfortunately observe a numerical underflow so that the determinant is outputted as zero (which it actually isn't).

How can I avoid this underflow?

One solution would be to obtain, instead of the determinant, the log of the determinant. However,

Any help is greatly appreciated.

Upvotes: 2

Views: 2424

Answers (2)

Alexander Shukaev
Alexander Shukaev

Reputation: 17021

There are 2 main options that come to my mind:

  1. The product of eigenvalues of square matrix is the determinant of this matrix, therefore a sum of logarithms of each eigenvalue is a logarithm of the determinant of this matrix. Assume det(A) = a and det(B) = b for compact notation. After applying aforementioned for 2 matrices A and B, we end up with log(a) and log(b), then actually the following is true:

    log(a + b) = log(a) + log(1 + e ^ (log(b) - log(a)))

    Yes, we get a logarithm of the sum. What would you do with it next? I don't know, depends on what you have to. If you have to remove logarithm by e ^ log(a + b) = a + b, then you might be lucky that the value of a + b does not underflow now, but in some cases it can still underflow as well.

  2. Perform clever preconditioning; there might be tons of options here, and you better read about them from some trusted sources as this is a serious topic. The simplest (and probably the cheapest ever) example of preconditioning for this particular problem could be to recall that det(c * A) = (c ^ n) * det(A), where A is n by n matrix, and to premultiply your matrix with some c, compute the determinant, and then to divide it by c ^ n to get the actual one.

Update

I thought about one more option. If on the last stages of #1 or #2 you still experience underflow too frequently, then it might be a good idea to increase precision specifically for these last operations, for example, by utilizing GNU MPFR.

Upvotes: 5

tmyklebu
tmyklebu

Reputation: 14205

You can use Householder elimination to get the QR decomposition of delta0. Then the determinant of the Q part is +/-1 (depending on whether you did an even or odd number of reflections) and the determinant of the R part is the product of the diagonal elements. Both of these are easy to compute without running into underflow hell---and you might not even care about the first.

Upvotes: 1

Related Questions