how to make a matrix invertible for sure, in R

Question

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?!

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.

This seems like a mathematics question, not an R or programming question. I suggest you ask this at [math.se] instead. — r2evans, Aug 23 '21 at 21:52
Can you add links to the posts you've read (with a short description of what you think each one is suggesting/why it's not working for you)? Have you tried `Matrix::nearPD()` and/or `sfsmisc::posdefify()`, `corpcor::make.positive.definite()`, `lqmm::make.positive.definite()`, `pracma::nearest_spd()` ... ? — Ben Bolker, Aug 23 '21 at 22:14
@BenBolker, I tried these all now. still same error. how come in R, a matrix with non zero determinant is not invertible? — Mathica, Aug 23 '21 at 22:40
Seems like it's time to produce the offending matrix as an [edit] to the question body. — IRTFM, Aug 23 '21 at 23:36
@BenBolker, You must definately be right! since excel invert the matrix easily. while r does not. I have a 45 by 45 matrix as a csv or rData file. I do not know how to insert the whole actual matrix here. and that matrix is result of a part of my huge code! — Mathica, Aug 24 '21 at 01:42
(I just found out Excel has native matrix inversion, I'm surprised!) please read the link I provided above, search for `dput()` ... — Ben Bolker, Aug 24 '21 at 01:49
@Ben Bolker, yeah, excel is getting more and more surprising! thanks for the dput link. however apparantly entering structure output from dput gives more than allowed 40000 character here. I am really confused about this. how a non zero determeninant matrix can not be inverted in R. — Mathica, Aug 24 '21 at 02:45
@BenBolker, thank you for your response. this is where I pasted my matrix : https://pastebin.com/TsUzuF4L — Mathica, Aug 24 '21 at 14:56

Ben Bolker · Accepted Answer · 2021-08-24T16:05:29.107

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

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

wonderful! so as far as I now realize, there s nothing to do with matrix itself ( like scaling, normalizing , ... ) to avoid this? — Mathica, Aug 24 '21 at 17:01
correct. All of these properties are (approximately) invariant to affine transformations (adding or multiplying by scalar (?) constants). The likelihood of a bad outcome (mistaken answer etc.) downstream depends entirely on the context of this computation and what you plan to do with the results. — Ben Bolker, Aug 24 '21 at 17:03

how to make a matrix invertible for sure, in R

1 Answers1