Numpy: Singular Matrix

user_9550385 0
Viewed 1040

I am trying to calculate the inverse of a matrix of the form X'X using NumPy as follows.

df = pd.read_csv('https://raw.githubusercontent.com/jianghaochu/data/master/x.csv')
X = np.array(df)
X.shape  # returns (92, 18)
np.linalg.det(np.matmul(X.transpose(), X))  # returns 8.029863818629298
np.linalg.matrix_rank(np.matmul(X.transpose(), X))  # returns 17

What confuses me is that X'X is not of full rank but the determinant is positive and NumPy can successfully calculate the inverse of X'X. However, if I take another matrix Y which is formed by a subset of the columns of X, then the determinant becomes zero and the inverse can not be calculated.

Y = X[:, [0, 12, 13, 14, 15, 16, 17]]
Y.shape  # returns (92, 7)
np.linalg.det(np.matmul(Y.transpose(), Y))  # returns 0.0
np.linalg.matrix_rank(np.matmul(Y.transpose(), Y))  # returns 16
np.linalg.inv(np.matmul(Y.transpose(), Y))  # numpy.linalg.LinAlgError: Singular matrix

numpy.linalg.LinAlgError: Singular matrix

In my opinion, the columns of Y are linear dependent. Hence, Y'Y is singular and its determinant is zero. Adding more columns to Y (like in X) should not lead to linear independence. Hence, I cannot understand how NumPy can calculate the inverse of X'X given Y'Y is singular.

I am calculating the OLS estimator using the data above. When I am using the columns (variables) in Y, I receive a singular error which I assume is due to perfect collinearity of the independent variables. When I add more independent variables in the model and get a matrix as X, I can successfully get the estimator using NumPy. It really confuses me since it suggests that adding more independent variables can solve perfect collinearity which really blows my head off.

I have been struggling with it for the entire afternoon and evening but still cannot get a clue. I also tried it in R and got the same result. I hope someone can guide me in some new directions. Any suggestions will be appreciated.

matrix-inverse linear-algebra numpy python
1 Answers

Check out the discussion here: numpy inverts a singular matrix.

Summarizing what has been said:

  • The reason you are getting such results is because numpy is using LU decomposition to calculate the inverse.
  • Such 'inverse' cannot be used to solve systems of linear equations.
  • This does not happen in Numpy 1.12.0
  • The condition number of your matrix is really high: np.linalg.cond(X.T@X) = 5.7294895077058016e+17

In addition to this, you can check eigenvalues - if one eigenvalue of the matrix is zero, its corresponding eigenvector is linearly dependent:

lambdas, V = np.linalg.eig(X.T@X)

lambdas[np.isclose(lambdas, 0)]
array([1.89103844e-15])
user_124115360
Related