x5g <- new("matrix.csr",
ra = c(300, 130, 5, 130, 330,
125, 10, 5, 125, 200, 70,
10, 70, 121.5, 1e30),
ja = c(1:3, 1:4, 1:4, 2:5),
ia = c(1L, 4L, 8L, 12L, 15L, 16L),
dimension = c(5L, 5L))
(m5g <- as.matrix(x5g)) # yes, is symmetric, and positive definite:
eigen(m5g, only.values=TRUE)$values # all positive (but close to singular)
ch5g <- chol(x5g)
str(ch5g) # --> the slots of the "matrix.csr.chol" class
mch5g <- as.matrix.csr(ch5g)
print.table(as.matrix(mch5g), zero.print=".") # indeed upper triagonal w/ positive diagonal
## x5 has even more extreme entry at [5,5]:
x5 <- x5g; x5[5,5] <- 2.9e32
m5 <- as.matrix(x5)
(c5 <- chol(m5))# still fine, w/ [5,5] entry = 1.7e16 and other diag.entries in (9.56, 17.32)
ch5 <- chol(x5) # --> warning "Replaced 3 tiny diagonal entries by 'Large'"
# gave error for a while
(mmc5 <- as.matrix(as.matrix.csr(ch5)))
# yes, these replacements were extreme, and the result is "strange'
## Solve the problem (here) specifying non-default singularity-tuning par 'tiny':
ch5. <- chol(x5, tiny = 1e-33)
(mmc5. <- as.matrix(as.matrix.csr(ch5.))) # looks much better.
## Indeed: R'R back-permuted *is* the original matrix x5, here m5:
(RtR <- crossprod(mmc5.)[ch5.@invp, ch5.@invp])
all.equal(m5, RtR, tolerance = 2^-52)
stopifnot(all.equal(m5, RtR, tolerance = 1e-14)) # on F38 Linux, only need tol = 1.25e-16
Run the code above in your browser using DataLab