5

I was attempting to calculate an OLS regression in R when I saw some strange things. The inverse of a square matrix does not exist if and only if the determinants is 0. Also, the matrix must be of full rank.

So not sure how the below is possible:

> dim(X)
[1] 20000    51
> det(t(X) %*% X)
[1] 3.863823e+161 #non-zero 

> solve(t(X) %*% X)
Error in solve.default(t(X) %*% X) :
system is computationally singular: reciprocal condition number = 3.18544e-17

Why is solve() throw an error when trying to calculate the inverse when we know the determinant is not zero? What am I missing here?

Checked that the matrix has full rank:

> qr(t(X) %*% X)$rank
[1] 51

But then just to test further I reassigned one of the X columns to the same value of another:

> X[,2] = X[,3]

Thus, two columns of the X matrix are now the same.

> qr(t(X) %*% X)$rank
[1] 50

We now can confirm the X'X matrix is not of full rank.

> det(t(X) %*% X)
[1] 1.634637e+138

But the determinant is still not equal to 0? How is this possible and what am I missing?

Rex.32
  • 51
  • 1
  • 5
  • 3
    As a general comment, you do not need to invert a matrix to solve an OLS regression. A much superior strategy is to explicitly solve the system of equations. – Matthew Drury Jan 13 '16 at 19:36
  • 3
    More elaboration on Matthew's point can be found here. http://stats.stackexchange.com/questions/160179/do-we-need-gradient-descent-to-find-the-coefficients-of-a-linear-regression-mode/164164#164164 But generally, a key lesson that you will learn when doing mathematics on computers is that formulas like the normal equation are sometimes of little *practical* utility, even if they're mathematically helpful. – Sycorax Jan 13 '16 at 19:38
  • @MatthewDrury@user777 I am aware that you can use SVD and the Moore-Penrose's pseudoinverse, or QR... this was not my point in asking the question. It was simply regarding the apparent violation of the invertible matrix theorem. – Rex.32 Jan 13 '16 at 20:36
  • 3
    @Roxxy.32 We know that, which is why we posted as comments instead of an answers. In general, the culture around here is to offer anything we think may be helpful, long or short term, in comments, while reserving the answers for addressing the precise question posed. No worries if we are telling you something you already know, because theres a good chance someone in the future will find your question, and find the comments helpful. – Matthew Drury Jan 13 '16 at 20:47
  • 2
    This may be helpful -- the numbers world from the point of view of a computer: https://ideas.repec.org/p/boc/scon14/22.html. A more dramatic explanation of the roundoff error is here: http://www.imdb.com/title/tt0151804/. – StasK Jan 13 '16 at 21:15

3 Answers3

9

My guess is that the numbers are too big (the determinant is large) and you're running into a computational problem.

I was able to replicate your error by running:

> X <- cbind(1,exp(rexp(100,rate=1/50)))
> det(t(X) %*% X)
[1] 5.156683e+126
> solve(t(X) %*% X)
> Error in solve.default...

The problem is numerical. You might be able to solve it by making some transformation of your $X$ matrix that makes the numbers smaller but allows you to work out what $\left(X'X\right)^{-1}$ is.

John
  • 415
  • 2
  • 6
  • seems to be the cas, thank you. It is counter intuitive since my X matrix is centered around 0 and has no values greater in abs value than 15, but seems the `t(X) %*% X` created some large number due to the sheer volume, however not larger than E^6. – Rex.32 Jan 13 '16 at 20:51
  • 1
    When the coefficients are on the order of $10^6$, the numbers involved in a $51\times 51$ determinant are on the order of $(10^6)^{51}=10^{306}$. If you wind up with a result of $10^{161}$, you may therefore have lost up to $145$ significant figures! – whuber Jan 14 '16 at 00:06
3

The method for determinant is different than the method for inverting a matrix. The determinant uses a lower upper decomposition. The determinant of a product is the product of determinants. The L is approximately very small and the U is approximately very large. At 16 point digit precision the very small number is rounded too large and the product explodes when it's actually 0.

I would trust the solve command. The matrix is singular. The r help says "you shouldn't use det for solving any problems".

AdamO
  • 52,330
  • 5
  • 104
  • 209
2

It looks like there's a similar question here, and I'd suggest a similar exploration.

What is the condition number of your matrix? Your matrix may be nearly singular, although I suspect that's unlikely.

What about the scale of $X$? What are its max values? Your determinant may be overflowing due to scaling issues, in which case you can decrease the values of the matrix by some constant factor.

I also agree with the commenters -- there's no need to explicitly invert a matrix to solve linear regression.

Community
  • 1
achompas
  • 181
  • 5