how to make a matrix invertible for sure, in R

808 Views Asked by At

this link has the dput output structure of my matrix.
https://pastebin.com/TsUzuF4L

Error in solve() : system is computationally singular: reciprocal condition number = 4.35295e-21 in R

I was wondering if there is any general method in R to make a matrix invertible for sure? any function?

I added the attribute tol=FALSE or tol = 1e-22 (compared to the number in error ) but I still get the same error.

ps. the reason I am bringing this here on stackexchange is, my matrix determinant is non zero ,but R gives the error above and believes my matrix is not invertible! how come?!

enter image description here

My matrix is 45 × 45. dput() output exceeds my limit of 40000 character on Stack Overflow, but to give an idea of what its figures are, I show part of it above.

1

There are 1 best solutions below

2
Ben Bolker On BEST ANSWER

tl;dr I was able to invert your matrix by setting tol=0 but it might not be a good idea.

When I get your matrix from the link you provided, I am able to work around the problems and invert the matrix, but I would suggest that you should be extremely cautious, as the warnings and errors are telling you that the computation is numerically unstable - you will probably get different answers on different operating systems, with different compilers, etc.. You can't trust these answers unless you have an independent way to verify them.

range(v <- eigen(M)$values)
## [1] -7.369628e-05  9.355280e+11
plot(abs(v), log="y", col = sign(v)+3, pch=16)

The matrix is not positive definite (which might be important depending on your application); the smallest eigenvalue is 16 orders of magnitude smaller than the largest ...

log-plot of eigenvalues

Matrix::rankMatrix(M)
## 19, with tolerance 1e-15
Matrix::rankMatrix(M, tol =0 )
## 45 with tolerance 0

When computed with the default tolerance, your matrix is reported as being rank-deficient, i.e. there are only 19 independent dimensions/columns (this corresponds to the number of eigenvalues above the big gap in the plot above)

We can compute the condition number:

Matrix::condest(M)
## $est: [1] 2.732966e+18

From Wikipedia:

As a rule of thumb, if the condition number κ(A) = 10^k, then you may lose up to k digits of accuracy on top of what would be lost to the numerical method due to loss of precision from arithmetic methods

In this case k=18 (draw your own conclusions) ...

When I compute the determinant, I get a very different value (but still non-zero).

det(M)
## [1] 2.568633e+35

I can invert the matrix if I tell R to ignore all of these warning signs by setting the tolerance to 0.

i <- solve(M, tol=0)

Depending on what you are doing, you might be interested in computing a pseudo-inverse that takes account of the (near) rank-deficiency of the matrix, e.g. using MASS::ginv().

Since the answers are likely to be highly dependent on system details, here is information from sessionInfo():

R Under development (unstable) (2021-07-30 r80684)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Pop!_OS 20.10

Matrix products: default
BLAS:   /usr/local/lib/R/lib/libRblas.so
LAPACK: /usr/local/lib/R/lib/libRlapack.so