6
$\begingroup$

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?

$\endgroup$
5
  • 3
    $\begingroup$ 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. $\endgroup$ Jan 13, 2016 at 19:36
  • 4
    $\begingroup$ 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. $\endgroup$
    – Sycorax
    Jan 13, 2016 at 19:38
  • $\begingroup$ @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. $\endgroup$
    – Rex.32
    Jan 13, 2016 at 20:36
  • 3
    $\begingroup$ @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. $\endgroup$ Jan 13, 2016 at 20:47
  • 2
    $\begingroup$ 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. $\endgroup$
    – StasK
    Jan 13, 2016 at 21:15

4 Answers 4

9
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – Rex.32
    Jan 13, 2016 at 20:51
  • 2
    $\begingroup$ 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! $\endgroup$
    – whuber
    Jan 14, 2016 at 0:06
3
$\begingroup$

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

$\endgroup$
2
  • 1
    $\begingroup$ The help page actually says "Often, computing the determinant is not what you should be doing to solve a given problem", which is not quite the same. The issue there is that a lot of statisticians try to compute determinants (and inverses) when they aren't needed. It's not saying there's a problem with the det algorithm (which is based on lapack's DGETRF) $\endgroup$ Apr 11, 2022 at 7:30
  • 1
    $\begingroup$ @ThomasLumley Good point. I suppose with multiple numerical algorithms, the answer to the question, "is my matrix approximately singular?" would be "it depends" in some cases. I suppose the point is: if you're fitting OLS with the default settings, what's most important is that solve can indeed solve. John's answer below points to the possibility of a numerical issue when the determinant is too big. $\endgroup$
    – AdamO
    Apr 11, 2022 at 15:24
2
$\begingroup$

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.

$\endgroup$
0
$\begingroup$

Ok, I think det is the one that's misleading here.

The "true" determinant is zero if the product of the eigenvalues of $X^TX$ is zero, which happens iff one of the individual eigenvalues is zero. Given computer arithmetic, the determinant will be computed as zero if one of the individual computed eigenvalues is exactly zero or if enough of them are very small that the computed product underflows. It takes a lot to underflow double precision, so we're talking really really small. .Machine$double.eps^20 doesn't underflow.

The matrix is truly uninvertible iff one of the individual eigenvalues is zero. Given computer arithmetic, the inverse will be detected as numerically singular if the estimated condition number, the ratio of the largest and smallest eigenvalues, is too large. The default threshold is the reciprocal of the condition number being smaller than machine epsilon, which is only $2^{-52}\approx 2\times 10^{-16}$. So it's a lot easier to get solve to give up on a matrix that to get det to underflow to zero.

@John's answer gives a matrix of rank 2 that has a non-zero computed determinant, because the non-zero eigenvalues are big and presumably the zero ones didn't exactly evaluate to zero. Your example isn't like that because it would have full rank at infinite precision, but it's presumably similar. The smallest eigenvalue is not zero, but it's less than machine epsilon times the largest eigenvalue.

As a final note, while solve and det just use Lapack, as all sensible people do, functions like lm and glm don't -- and they have a much stricter tolerance for singular matrices, because typically a double-precision design matrix that someone hasn't deliberately set up as a numerical analysis exercise is either actually singular or has a reciprocal condition number much larger than machine epsilon. And if it does fall in the gap, the user probably needs to know. The tolerance (in qr(,LAPACK=FALSE)) is $10^{-7}$. So, the numerical rank as computed by qr can be zero when solve still works, and that's deliberate and for good reasons. (I mean, on top of the fact that you're probably using qr on $X$ rather than $X^TX$)

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.