matrix.csr.chol-class.RdA class of objects returned from Ng and Peyton's (1993) block sparse Cholesky algorithm.
Objects may be created by calls of the form new("matrix.csr.chol",
...), but typically result from
chol(<matrix.csr>).
nrow:an integer, the number of rows of the
original matrix, or in the linear system of equations.
nnzlindx:Object of class numeric, number of non-zero elements in lindx
nsuper:an integer, the number of supernodes of the decomposition
lindx:Object of class integer, vector of integer
containing, in column major order, the row subscripts of the non-zero entries
in the Cholesky factor in a compressed storage format
xlindx:Object of class integer, vector of integer of pointers for lindx
nnzl:of class "numeric", an integer, the number of non-zero
entries, including the diagonal entries, of the Cholesky factor stored in lnz
lnz:a numeric vector of the entries of the Cholesky factor
xlnz:an integer vector, the column pointers for the Cholesky factor stored in lnz
invp:inverse permutation vector, integer
perm:permutation vector, integer
xsuper:Object of class integer, array containing the supernode partioning
det:numeric, the determinant of the Cholesky factor
log.det:numeric, the log determinant of the Cholesky factor
ierr:an integer, the error flag (from Fortran's src/chol.f)
time:numeric, unused (always 0.) currently.
Note that the perm and notably invp maybe important to back
permute rows and columns of the decompositions, see the Examples, and our
chol help page.
signature(x = "matrix.csr.chol",
upper.tri=TRUE): to get the sparse ("matrix.csr") upper
triangular matrix corresponding to the Cholesky decomposition.
signature(r = "matrix.csr.chol"): for computing
\(R^{-1} b\) when the Cholesky decomposition is \(A = R'R\).
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:
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 300 130 5 0.0 0e+00
#> [2,] 130 330 125 10.0 0e+00
#> [3,] 5 125 200 70.0 0e+00
#> [4,] 0 10 70 121.5 0e+00
#> [5,] 0 0 0 0.0 1e+30
eigen(m5g, only.values=TRUE)$values # all positive (but close to singular)
#> [1] 1.000000e+30 4.846141e+02 2.642875e+02 1.433734e+02 5.922500e+01
ch5g <- chol(x5g)
str(ch5g) # --> the slots of the "matrix.csr.chol" class
#> Formal class 'matrix.csr.chol' [package "SparseM"] with 15 slots
#> ..@ nrow : int 5
#> ..@ nnzlindx: int 7
#> ..@ nsuper : int 3
#> ..@ lindx : int [1:7] 1 2 4 5 3 4 5
#> ..@ xlindx : int [1:6] 1 2 5 8 11 11
#> ..@ nnzl : int 10
#> ..@ lnz : num [1:10] 1.00e+15 1.10e+01 6.35 9.07e-01 1.73e+01 ...
#> ..@ xlnz : int [1:6] 1 2 5 8 10 11
#> ..@ invp : int [1:5] 3 5 4 2 1
#> ..@ perm : int [1:5] 5 4 1 3 2
#> ..@ xsuper : int [1:6] 1 2 3 6 1 1
#> ..@ det : num 3.3e+19
#> ..@ log.det : num 44.9
#> ..@ ierr : int 0
#> ..@ time : num 0
mch5g <- as.matrix.csr(ch5g)
print.table(as.matrix(mch5g), zero.print=".") # indeed upper triagonal w/ positive diagonal
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 1.000000e+15 . . . .
#> [2,] . 1.102270e+01 . 6.350529e+00 9.072184e-01
#> [3,] . . 1.732051e+01 2.886751e-01 7.505553e+00
#> [4,] . . . 1.263279e+01 9.267311e+00
#> [5,] . . . . 1.367335e+01
## 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)
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 17.32051 7.505553 0.2886751 0.0000000 0.000000e+00
#> [2,] 0.00000 16.542874 7.4251509 0.6044899 0.000000e+00
#> [3,] 0.00000 0.000000 12.0326140 5.4445003 0.000000e+00
#> [4,] 0.00000 0.000000 0.0000000 9.5651455 0.000000e+00
#> [5,] 0.00000 0.000000 0.0000000 0.0000000 1.702939e+16
ch5 <- chol(x5) # --> warning "Replaced 3 tiny diagonal entries by 'Large'"
#> Warning: In chol(<matrix.csr>, *): Replaced 3 tiny diagonal entries by 'Large'
# gave error for a while
(mmc5 <- as.matrix(as.matrix.csr(ch5)))
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 1.702939e+16 0e+00 0.00000 0.000000e+00 0.000000e+00
#> [2,] 0.000000e+00 1e+64 0.00000 7.000000e-63 1.000000e-63
#> [3,] 0.000000e+00 0e+00 17.32051 2.886751e-01 7.505553e+00
#> [4,] 0.000000e+00 0e+00 0.00000 1.000000e+64 1.228333e-62
#> [5,] 0.000000e+00 0e+00 0.00000 0.000000e+00 1.000000e+64
# 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.
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 1.702939e+16 0.0000 0.00000 0.0000000 0.0000000
#> [2,] 0.000000e+00 11.0227 0.00000 6.3505290 0.9072184
#> [3,] 0.000000e+00 0.0000 17.32051 0.2886751 7.5055535
#> [4,] 0.000000e+00 0.0000 0.00000 12.6327926 9.2673109
#> [5,] 0.000000e+00 0.0000 0.00000 0.0000000 13.6733526
## Indeed: R'R back-permuted *is* the original matrix x5, here m5:
(RtR <- crossprod(mmc5.)[ch5.@invp, ch5.@invp])
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 300 130 5 0.0 0.0e+00
#> [2,] 130 330 125 10.0 0.0e+00
#> [3,] 5 125 200 70.0 0.0e+00
#> [4,] 0 10 70 121.5 0.0e+00
#> [5,] 0 0 0 0.0 2.9e+32
all.equal(m5, RtR, tolerance = 2^-52)
#> [1] TRUE
stopifnot(all.equal(m5, RtR, tolerance = 1e-14)) # on F38 Linux, only need tol = 1.25e-16