4

I am trying to estimate the mean value of log(det(AAT)+1) in Python. My simple code works fine until I get to 17×17 matrices at which point it gives me a math error. Here is the code:

iter = 10000
for n in xrange(1,20):
    h = n
    dets = []
    for _ in xrange(iter):
        A = (np.random.randint(2, size=(h,n)))*2-1
        detA_Atranspose = np.linalg.det(np.dot(A, A.transpose()))
        try:
            logdetA_Atranspose = math.log(detA_Atranspose+1,2)
        except ValueError:
            print "Ooops!", n,detA_Atranspose
        dets.append(logdetA_Atranspose)
    print np.mean(dets)

A is supposed to be a matrix with elements that are either -1 or 1.

What am I doing wrong and how can it be fixed? What is special about 17?

Nayuki
  • 17,911
  • 6
  • 53
  • 80
Simd
  • 19,447
  • 42
  • 136
  • 271

1 Answers1

2

For the formula in the title ( previously logdet(AA^T) ):

det(AA^T) for some random As can simply be 0. The function would then fail because computing log(0) is invalid.

Note that in theory det(AA^T) cannot be negative, since AA^T is a positive semi-definite matrix (which means that all eigenvalues are non-negative and implies that det >= 0).

For the formula in the code ( logdet(1+AA^T) )

You should probably use numpy.linalg.slogdet() and compute slogdet(1+A.dot(A.T))

From its documentation:

"Compute the sign and (natural) logarithm of the determinant of an array.

If an array has a very small or very large determinant, then a call to det may overflow or underflow. This routine is more robust against such issues, because it computes the logarithm of the determinant rather than the determinant itself."

Tomer Levinboim
  • 992
  • 12
  • 18
  • You're right about the +1 in log(1+AA^T). I guess I took it from the title... so, it can only happen due to overflow. – Tomer Levinboim Mar 06 '16 at 11:11
  • 1
    The proposed use of `slogdet` should convert the computation to floating point, which also will repair the integer overflow problem. However, with a rather small probability, the determinant can still be zero. – Lutz Lehmann Mar 06 '16 at 11:20
  • "det(AA^T) cannot be negative"... it does appear to be though according to the code. – Simd Mar 06 '16 at 13:45
  • @eleanora only due to overflow, presumably. – Cimbali Mar 06 '16 at 19:44
  • slogdet(1+A.dot(A.T)) doesn't give 0 when A is singular sadly so I get a lot of `-inf` values. – Simd Mar 06 '16 at 20:00
  • @Cimbali: No, only many medium sized eigenvalues that amplify floating point noise at zero eigenvalues. Like degenerate noisy eigenvalue times remaining product of eigenvalues equals `-1e-15*16^14=-2^-50*2^56=-2^6`, as can happen in the case `n=16`. – Lutz Lehmann Mar 06 '16 at 20:15
  • @eleanora:Just cut off the negative determinants at zero as obvious results of a noisy singular matrix. – Lutz Lehmann Mar 06 '16 at 20:16
  • 1
    This could be another solution to your original problem: compute logdet(1 + AA^T + pI), where p is a small scalar (say 1e-6), and I is the identity matrix. This should get rid of instability issues at the expense of accuracy. – Tomer Levinboim Mar 06 '16 at 20:52
  • @LutzL are you sure? Because on a 64 bit signed integer overflow would occur at 2^63, and `ln(2^63) = 43.6`, not impossibly far off values already observed... I don't know for a fact whether the internals of `numpy.linalg.det` are done in floating point or integer math, but the documentation does say it is subject to overflows. – Cimbali Mar 07 '16 at 12:15
  • Though I guess it doesn't make much sense to use integer math if you may expect floating point values in your input matrices. – Cimbali Mar 07 '16 at 12:16
  • @Cimbali: Yes, that is very suggestive, but the problem stays the same if you force floating point values via `A = (np.random.randint(2, size=(h,n)))*2-1.0`. Also, the wrongly negative values are not negative integers, but have non-trivial fractional parts. – Lutz Lehmann Mar 07 '16 at 17:21