[R] impossible to invert a spam-object, but possible when it's a matrix-object

Berend Hasselman bhh at xs4all.nl
Mon Feb 11 16:01:33 CET 2013


On 11-02-2013, at 15:00, "Camarda, Carlo Giovanni" <Camarda at demogr.mpg.de> wrote:

> Dear Uwe,
> 
> thanks for the response.
> I knew my matrix was almost singular, but is there any way to implement the algorithm is solve(base) with a spam-matrix?
> 
> As workaround, I might add a ridge penalty:
> 
>    ApP <- A+10^-1*diag.spam(nrow(A))
>    solve(ApP)
> 
> but this would completely jeopardize other parts of my model.


I've inspected your  matrix with package Matrix.
It doesn't seem catastrophically ill-conditioned.
It's not positive definite.

It seems that spam is issuing a misleading error message.

Try this

ent <- c(2312.12324929972,-2000,1000,-2000,1000,-2000,
         5031.91011235955,-2000,-2000,1000,-1,1000,-2000,
         2049.8595505618,-2000,1000,-1,-2000,5036.89635854342,
         -2000,1000,-2000,1,-2000,-2000,8119.66292134831,-2000,
         -2000,-2000,1000,-2000,5058.83426966292,-2000,-1,1000,
         -2000,2051.85434173669,-2000,1000,1,1000,-2000,-2000,
         5043.87640449438,-2000,1,1000,-2000,1000,-2000,
         2110.68820224719,-1,1,0,-1,1,-1,1)
rr <- as.integer(c(1,6,12,18,24,29,35,41,47,52,55,57,59))
cc <- as.integer(c(1,2,3,4,7,1,2,3,5,8,10,1,2,3,6,9,11,1,
                   4,5,6,7,10,2,4,5,6,8,3,4,5,6,9,12,1,4,
                   7,8,9,11,2,5,7,8,9,12,3,6,7,8,9,2,4,10,3,7,6,8))
di <- as.integer(c(12,12))

library(Matrix)
A <- sparseMatrix(x=ent,p=rr-1,j=cc)
A
condest(A)
det(A)
solve(A)
kappa(A)

# not pd
# chol(A)
isSymmetric(A)
# chol(A) # ==> not positive definite
chol(A+diag(1e-1,nrow(A))) # this works but isn't what you want

A.m <- as.matrix(A)
chol(A.m,pivot=TRUE) # without pivot: not positive definite
# has rank 9 and warning message

B <- as.matrix(A) 
condest(B)
rcond(B)
det(A)
solve(B)

all.equal(as.matrix(A %*% solve(B)), diag(12), check.attributes=FALSE)
all.equal(as.matrix(B %*% solve(A)), diag(12), check.attributes=FALSE)

Berend



More information about the R-help mailing list