matrix.csr.chol-class function

Class "matrix.csr.chol" (Block Sparse Cholesky Decomposition)

A class of objects returned from Ng and Peyton's (1993) block sparse Cholesky algorithm. class

Objects from the Class

Objects may be created by calls of the form new("matrix.csr.chol", ...), but typically result from chol(<matrix.csr>).

Slots

• 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.

Details

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.

Methods

• as.matrix.csr: signature(x = "matrix.csr.chol", upper.tri=TRUE): to get the sparse ("matrix.csr") upper triangular matrix corresponding to the Cholesky decomposition.
• backsolve: signature(r = "matrix.csr.chol"): for computing $R^{-1} b$ when the Cholesky decomposition is $A = R'R$.

See Also

Base 's chol and SparseM's chol, notably for examples; backsolve

Examples

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