# Matrix inverse not able to be calculated while determinant is non-zero

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)
 20000    51
> det(t(X) %*% X)
 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  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
 50


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

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


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

• 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
• More elaboration on Matthew's point can be found here. stats.stackexchange.com/questions/160179/… 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. – Reinstate Monica 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
• This may be helpful -- the numbers world from the point of view of a computer: ideas.repec.org/p/boc/scon14/22.html. A more dramatic explanation of the roundoff error is here: imdb.com/title/tt0151804. – StasK Jan 13 '16 at 21:15

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

• 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
• 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 0:06

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

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.