Search code examples
rmatrixmatrix-inversesingular

how to make a matrix invertible for sure, in R

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.

Solution

  • 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