[R] Can't invert matrix

[Bill.Venables at csiro.au]

What you show below is only a representation of the matrix to 7dp.  If you look at that, though, the condition number is suspiciously large (i.e. the matrix is very ill-conditioned):

> txt <- textConnection("
+   0.99252358 0.93715047 0.7540535 0.4579895
+   0.01607797 0.09616267 0.2452471 0.3088614
+   0.09772828 0.58451468 1.4907090 1.8773815
+  -0.01000000 0.00000000 0.0900000 0.1700000
+ ")
> mat <- matrix(scan(txt, quiet = TRUE), ncol = 4, byrow = TRUE)
> close(txt)
> 
> with(svd(mat), d[1]/d[4])
[1] 82473793

With this representation, though, the matrix does seem to allow inversion on a windows 32 bit machine:

> solve(mat)
           [,1]      [,2]     [,3]        [,4]
[1,]  1.4081192  -7826819  1287643  21.5115344
[2,] -0.3987761  18796863 -3092403 -25.7757002
[3,] -0.1208272 -18836093  3098860  -0.9159397
[4,]  0.1467979   9511648 -1564829   7.6326466
> 

and horrible as it is (look closely at columns 2 and 3) it does check out pretty well:

> mat %*% solve(mat)
              [,1]          [,2]          [,3]          [,4]
[1,]  1.000000e+00 -3.060450e-10  1.108447e-10 -1.025872e-15
[2,]  3.571091e-18  1.000000e+00 -5.269385e-11 -1.840975e-16
[3,]  8.201989e-17  2.559318e-09  1.000000e+00 -1.284563e-15
[4,] -3.203479e-18 -8.867573e-11  1.813305e-11  1.000000e+00
> 

but a generalized inverse (with default tolerances) is very different

> library(MASS)
> ginv(mat)
            [,1]       [,2]       [,3]      [,4]
[1,]  1.27299552 -0.2800302 -1.7022000  16.38665
[2,] -0.07426349  0.1804989  1.0971974 -13.46780
[3,] -0.44601712  0.2273789  1.3821521 -13.24953
[4,]  0.31100880 -0.1368457 -0.8318576  13.86073
> 

which emphasises the delicacy of dealing with an ill-conditioned matrix.  Incidently this also checks out fairly well according to the definition of a ginverse:

> range(mat - mat %*% ginv(mat) %*% mat)
[1] -2.132894e-08  2.128452e-08
>

When dealing with numerical matrices you have to be prepared for the unexpected.

Bill Venables.
 

-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf Of Noah Silverman
Sent: Sunday, 21 November 2010 2:56 PM
To: r-help at r-project.org
Subject: [R] Can't invert matrix

Hi,

I'm trying to use the solve() function in R to invert a matrix.  I get
the following error, "Lapack routine dgesv: system is exactly singular"

However, My matrix doesn't appear to be singular. 

            [,1]       [,2]      [,3]      [,4]
[1,]  0.99252358 0.93715047 0.7540535 0.4579895
[2,]  0.01607797 0.09616267 0.2452471 0.3088614
[3,]  0.09772828 0.58451468 1.4907090 1.8773815
[4,] -0.01000000 0.00000000 0.0900000 0.1700000


Can anyone help me understand what is happening here?

Thanks!

-N

______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.