Search code examples
pythonnumpymatrixmatrix-inversedeterminants

numpy: Possible for zero determinant matrix to be inverted?


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

Solution

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