Avoid numerical underflow when obtaining determinant of large matrix in Eigen

Question

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,

I do not know how to do this;
how could I then add up these logs?

Any help is greatly appreciated.

You could multiply the matrix by a scalar before calculating the determinant, say 1000 * A. Then divide the result by 1000. Or you could work with bigger variables (long double). — lucas92, Dec 11 '13 at 21:11
@lucas92 1000 would be a poor choice, as multiplication by 1000 is generally inexact. An appropriate power of two instead would avoid making the computation less precise, even when it is unnecessary, as long as no overflow occurs. — Pascal Cuoq, Dec 11 '13 at 21:21

Alexander Shukaev · Accepted Answer · 2013-12-12T09:22:17.510

5

There are 2 main options that come to my mind:

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.
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.

edited Dec 12 '13 at 09:22

answered Dec 11 '13 at 21:16

Alexander Shukaev

16,674
8
70
85

The first solution sounds pretty appealing. Do you also know how I then can get the sum of the values (given that I only have their logs)? – Henrik Dec 11 '13 at 21:21
Do you mean that you have, for instance, matrices `A` and `B` and their corresponding logarithms of determinants, i.e. `log(det(A))` and `log(det(B))`, and now you actually want to know `det(A) + det(B)`? If that's what you want, then I don't see any problem with extracting, for example, `det(A)` from `log(det(A))` by doing `e ^ log(det(A))`... Otherwise, I misunderstood your intent with the sum. – Alexander Shukaev Dec 11 '13 at 21:28
Yes, that is what I intent. But the problem still is that `e ^ log(det(A))` is 0 due to numerical underflow so I have to avoid it. However, I think the solution with preconditioning is actually even easier. As I am no computer scientist nor mathematician, what would be a "trusted source" to check the formula you gave? Any book on matrix algebra? – Henrik Dec 11 '13 at 21:32
That's the problem, I can't remember the book. But it would be good to start with basic terminology from **[Wikipedia](http://en.wikipedia.org/wiki/Preconditioner)** as there are many different preconditioners which serve different purposes, and are mostly targeted at solution of linear systems. For example, Incomplete Cholesky Decomposition is definitely not for your case. In fact, I think you should be fine with the simplest method that I've already outlined for you in the answer. Better try it first as digging into preconditioning will eat quite some time from you (it's really nasty topic). – Alexander Shukaev Dec 11 '13 at 21:44
1

Updated, take a look. – Alexander Shukaev Dec 11 '13 at 21:55
Awesome, thanks a lot. I am willing to accept, but I will wait until I have this successfully implemented (could take some time). – Henrik Dec 11 '13 at 21:59

score 1 · Answer 2 · answered Dec 12 '13 at 03:58

1

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.

answered Dec 12 '13 at 03:58

tmyklebu

13,915
3
28
57

There is no guarantee that the determinant of `R` will not underflow. It heavily depends on the underlying problem (matrix). – Alexander Shukaev Dec 12 '13 at 09:06
1

@Haroogan: Correct. But it becomes trivial to find the logarithm of the determinant of R; it's the sum of the logarithms of the diagonal entries. Or you can use a similar trick to the one you use for geometric means to get the significand in a `double` and the exponent in a large enough integral type. – tmyklebu Dec 12 '13 at 09:10

Avoid numerical underflow when obtaining determinant of large matrix in Eigen

2 Answers2

Update