2

I have a estimate of covariance matrix

I want to take the inverse of this matrix, R gives me the following error

A=[ 3529861.470  8785861.47  6920.344 17120.34;
     8785861.470 26209861.47 17120.344 51920.34;
    6920.344    17120.34    14.000    34.00;
    17120.344    51920.34    34.000   104.00]

"Error in solve.default(l) : system is computationally singular: reciprocal condition number = 2.14511e-22".

however Matlab does the inverse calculation without giving any error message. Does anyone know the reason why R is giving me error? Is the any other way to calculate inverse

Spandyie
  • 914
  • 2
  • 11
  • 23
  • 1
    Matlab gives a condition number (`cond(A)`) of 9.5e12, so the matrix is actually close to singular. Perhaps R just has a stricter threshold than Matlab to declare a matrix as close to singular – Luis Mendo Feb 13 '15 at 22:48
  • reciprocal condition number that R gives me is 2.14511e-22 . isnt there a way to get around this situation in R? – Spandyie Feb 13 '15 at 22:50
  • I don't know R, but in linear algebra if you have large conditional number, algorithm for inverting should be a little bit different (I mean, you cannot calculate it with determinants, since you will get large numbers, you should use something like gauss elimination). May be you should search for different options of inverse matrix calculation in R – Mikhail Genkin Feb 13 '15 at 22:53
  • Get around what? Your matrix _is_ close to singular, so whatever you compute is subject to potentially large errors. Perhaps R has some setting to relax its close-to-singularness criterion. Are you using double-precision floats? – Luis Mendo Feb 13 '15 at 22:53
  • 1
    Your example isn't reproducible. `solve()` has no trouble inverting this matrix. – Richard Border Feb 13 '15 at 22:56
  • 1
    That's because solve use gauss elimination avoiding determinant calculation – Mikhail Genkin Feb 13 '15 at 22:56
  • 1
    @Spandy MATLAB's `rcond` returns `7.5E-14`. Are you sure that MATLAB and R have the same matrix? – TroyHaskin Feb 13 '15 at 22:58
  • The reason I am confused because the determinant of matrix isnt even very small. I meant isnt there a way to relax the thresholding in R? @LuisMendo – Spandyie Feb 13 '15 at 22:59
  • @TroyHaskin I get the determinant of -6.7604e+07 and Inverse =1.0e+05 * [ -0.0000 0.0000 0.0054 -0.0035; 0.0000 -0.0000 -0.0033 0.0021; 0.0054 -0.0033 -2.2911 1.4946; -0.0035 0.0021 1.4945 -0.9749] – Spandyie Feb 13 '15 at 23:01
  • The determinant can be very large and still the matrix is ill-conditioned (close to singular). For examply, take any matrix that is close to singular. Now multiply it by 10^6. There you go, the determinant is 10^6 larger and your matrix is as close to singular as before. To see if it's close to singular use the condition number, not the determinant. – Luis Mendo Feb 13 '15 at 23:01
  • 1
    @Spandy The determinant is [not suitable as a measure of conditioning](http://scicomp.stackexchange.com/a/1330). – TroyHaskin Feb 13 '15 at 23:02
  • 2
    `solve` inverts this matrix with no errors or warnings for me in R as well. – joran Feb 13 '15 at 23:16
  • can you give us `dput(A)` so we can see exactly what's there? – Ben Bolker Feb 14 '15 at 00:55

1 Answers1

5
A <- matrix(c(3529861.470,8785861.47,6920.344,17120.34,
       8785861.470,26209861.47,17120.344,51920.34,
   6920.344,17120.34,14.000,34.00,
   17120.344,51920.34,34.000,104.00),
     nrow=4,byrow=TRUE)

solve(A) ## works on my system

##              [,1]         [,2]         [,3]        [,4]
## [1,]   -1.2515442    0.7617239     535.4871   -349.3141
## [2,]    0.7617072   -0.4635922    -325.9051    212.5957
## [3,]  535.4884664 -325.9130639 -229114.3516 149458.2734
## [4,] -349.3061387  212.5955306  149454.4973 -97492.7335

eigen(A)$values
## [1]  2.921525e+07  5.245875e+05  1.440703e+00 -3.061760e-06

rcond(A)  ## condition number
## [1] 7.516179e-14

You should be able to adjust the tol parameter if you are having trouble inverting a matrix, but this is of course done at your own risk -- you're overriding a warning that the matrix operation may be numerically unstable.

This is with

R Under development (unstable) (2015-02-11 r67792)
Platform: i686-pc-linux-gnu (32-bit)
Running under: Ubuntu precise (12.04.5 LTS)
Ben Bolker
  • 211,554
  • 25
  • 370
  • 453