CODEWITHSUNDEEP
rmatrixinversematrix-inverse
Mayank Raj
Mayank Raj

Reputation: 121

Why does solve in R give weird inverse?

M matrix

     1000          793504755
793504755    629649829362825

This is my matrix M in R solve(M) gives this error : 

solve(M)
Error in solve.default(M) : 
  system is computationally singular: reciprocal condition number = 8.36282e-20

while MMINV in excel gives the desired output for this 2 x 2 matrix

Is R computing inverse using some decomposition algorithm? Why is it not giving output for my matrix? The matrix determinant is also not zero. in fact it is a large value.

Upvotes: 1

Views: 2767

Answers (2)

James
James

Reputation: 66834

You are having trouble because your matrix is ill-conditioned. You can solve it using the Matrix package without having to worry about intermediate steps:

library(Matrix)
solve(Matrix(M))
2 x 2 Matrix of class "dsyMatrix"
              [,1]          [,2]
[1,]  1.899097e+04 -2.393303e-02
[2,] -2.393303e-02  3.016117e-08

But there is some numerical imprecision in the result:

M %*% solve(Matrix(M))
2 x 2 Matrix of class "dgeMatrix"
            [,1] [,2]
[1,] 1.000000004    0
[2,] 0.001953125    1

This is because of the reciprocal condition number is smaller than the numerical precision afforded by doubles.

rcond(M)
[1] 8.362817e-20
storage.mode(M)
[1] "double"
.Machine$double.eps
[1] 2.220446e-16

However, the package gmp offers exact precision using rational numbers to give you the result you want:

library(gmp)
solve(as.bigz(M))
Big Rational ('bigq') 2 x 2 matrix:
     [,1]                  [,2]              
[1,] 932814562019/49118837 -5877813/245594185
[2,] -5877813/245594185    40/1326208599     
M %*% solve(matrix(as.bigz(c(M)),2))
Big Rational ('bigq') 2 x 2 matrix:
     [,1] [,2]
[1,] 1    0   
[2,] 0    1

Upvotes: 4

Roland
Roland

Reputation: 132706

This suffers from numerical difficulties. You can either use the "manual" method suggested by David or use the Choleski decomposition to get a somewhat reasonable estimate:

chol2inv(chol(M))
#              [,1]          [,2]
#[1,]  1.899097e+04 -2.393303e-02
#[2,] -2.393303e-02  3.016117e-08

chol2inv(chol(M)) %*% M
#             [,1]        [,2]
#[1,] 1.000000e+00 0.001953125
#[2,] 3.552714e-15 1.000000000

However, I'd suggest to avoid calculating the inverse of this matrix.

Upvotes: 3

Related Questions