3

By definition, a square matrix that has a zero determinant should not be invertible. However, for some reason, after generating a covariance matrix, I take the inverse of it successfully, but taking the determinant of the covariance matrix ends up with an output of 0.0.

What could be potentially going wrong? Should I not trust the determinant output, or should I not trust the inverse covariance matrix? Or both?

Snippet of my code:

cov_matrix = np.cov(data)
adjusted_cov = cov_matrix + weight*np.identity(cov_matrix.shape[0]) # add small weight to ensure cov_matrix is non-singular
inv_cov = np.linalg.inv(adjusted_cov) # runs with no error, outputs a matrix
det = np.linalg.det(adjusted_cov) # ends up being 0.0
kk415kk
  • 1,227
  • 1
  • 14
  • 30
  • Could you show us an example input and output? Your comment implies that you are attempting to make the matrix non-singular, which would mean that the matrix should have a valid inverse. – BobChao87 Feb 25 '15 at 06:40
  • Hmm, my dataset consists of 28x28 images where I use raw pixels as features, so the input output would be too large to copy, I think (784x487 covariance matrices). I believe adding the weight does result in having a valid inverse - however, why is the determinant zero then? If I don't add the weight, I get the singular lin algebra error. Would the rounding error apply to the determinant, where it's rounding to 0.0? – kk415kk Feb 25 '15 at 06:45
  • When you write `nv_cov = np.linalg(adjusted_cov)` do you mean `nv_cov = np.linalg.inv(adjusted_cov)` ? and should `det = np.linalg(adjusted_cov)` really be `det = np.linalg.det(adjusted_cov)` ? – gboffi Feb 25 '15 at 07:13
  • @kk415kk Since the numbers of the covarient matrix are small (according to your comment), it's possible you are running into an issue with not rounding, but rather precision. The product of 784 numbers with an average size of even just 0.3 is too small to be held in a Python float. If numpy has a method of computing the inverse of your matrix without computing the determinant, then it could possibly lead to the behaviour you describe. – BobChao87 Feb 25 '15 at 07:45
  • Ahhh okay. The context of my problem is that I'm attempting to model these images (sorted into classes) as a multivariate Gaussian distribution, which requires me to calculate the covariance matrix and its inverse. Is there a better way to do this if precision is an issue? – kk415kk Feb 25 '15 at 08:02
  • @kk415kk It is possible that numpy is modeling this correcting. It may only be the determinant that is failing. To test this, I recommend seeing if the product of `adjusted_cov` and `inv_cov` are (approximately) the identity matrix. If that's the case, then numpy probably isn't far off, especially if happens for multiple tests cases. – BobChao87 Feb 25 '15 at 08:09
  • It looks like the cov matrix multiplied to the inverse is resulting in the identity matrix. Is there a way to more accurately represent the determinant for my scenario then? – kk415kk Feb 25 '15 at 08:50
  • If you require the value of the determinant for other reasons, I'd suggest seeing if you can integrate numpy with the `decimal.Decimal` class which supports arbitrary precision. It will add a _lot_ of overhead to any computations done with that method, so I'd recommend only using `Decimal` when absolutely necessary. – BobChao87 Feb 25 '15 at 17:07

1 Answers1

3

The numerical inversion of matrices does not involve computing the determinant. (Cramer's formula for the inverse is not practical for large matrices.) So, the fact that determinant evaluates to 0 (due to insufficient precision of floats) is not an obstacle for the matrix inversion routine.

Following up on the comments by BobChao87, here is a simplified test case (Python 3.4 console, numpy imported as np)

A = 0.2*np.identity(500)
np.linalg.inv(A)

Output: a matrix with 5 on the main diagonal, which is the correct inverse of A. 

np.linalg.det(A)

Output: 0.0, because the determinant (0.2^500) is too small to be represented in double precision.

A possible solution is a kind of pre-conditioning (here, just rescaling): before computing the determinant, multiply the matrix by a factor that will make its entries closer to 1 on average. In my example, np.linalg.det(5*A) returns 1.

Of course, using the factor of 5 here is cheating, but np.linalg.det(3*A) also returns a nonzero value (about 1.19e-111). If you try np.linalg.det(2**k*A) for k running through modest positive integers, you will likely hit one that will return nonzero. Then you will know that the determinant of the original matrix was approximately 2**(-k*n) times the output, where n is the matrix size.