Cholesky {Matrix} R Documentation

## Cholesky Decomposition of a Sparse Matrix

### Description

Computes the Cholesky (aka “Choleski”) decomposition of a sparse, symmetric, positive-definite matrix. However, typically chol() should rather be used unless you are interested in the different kinds of sparse Cholesky decompositions.

### Usage

Cholesky(A, perm = TRUE, LDL = !super, super = FALSE, Imult = 0, ...)


### Arguments

 A sparse symmetric matrix. No missing values or IEEE special values are allowed. perm logical scalar indicating if a fill-reducing permutation should be computed and applied to the rows and columns of A. Default is TRUE.
 LDL logical scalar indicating if the decomposition should be computed as LDL' where L is a unit lower triangular matrix. The alternative is LL' where L is lower triangular with arbitrary diagonal elements. Default is TRUE. Setting it to NA leaves the choice to a CHOLMOD-internal heuristic. super logical scalar indicating if a supernodal decomposition should be created. The alternative is a simplicial decomposition. Default is FALSE. Setting it to NA leaves the choice to a CHOLMOD-internal heuristic. Imult numeric scalar which defaults to zero. The matrix that is decomposed is A+m*I where m is the value of Imult and I is the identity matrix of order ncol(A). ... further arguments passed to or from other methods.

### Details

This is a generic function with special methods for different types of matrices. Use showMethods("Cholesky") to list all the methods for the Cholesky generic.

The method for class dsCMatrix of sparse matrices — the only one available currently — is based on functions from the CHOLMOD library.

Again: If you just want the Cholesky decomposition of a matrix in a straightforward way, you should probably rather use chol(.).

Note that if perm=TRUE (default), the decomposition is

A = P' \tilde{L} D \tilde{L}' P = P' L L' P,

where L can be extracted by as(*, "Matrix"), P by as(*, "pMatrix") and both by expand(*), see the class CHMfactor documentation.

Note that consequently, you cannot easily get the “traditional” cholesky factor R, from this decomposition, as

R'R = A = P'LL'P = P'\tilde{R}'\tilde{R} P = (\tilde{R}P)' (\tilde{R}P),

but \tilde{R}P is not triangular even though \tilde{R} is.

### Value

an object inheriting from either "CHMsuper", or "CHMsimpl", depending on the super argument; both classes extend "CHMfactor" which extends "MatrixFactorization".

In other words, the result of Cholesky() is not a matrix, and if you want one, you should probably rather use chol(), see Details.

### References

Yanqing Chen, Timothy A. Davis, William W. Hager, and Sivasankaran Rajamanickam (2008) Algorithm 887: CHOLMOD, Supernodal Sparse Cholesky Factorization and Update/Downdate. ACM Trans. Math. Softw. 35, 3, Article 22, 14 pages. doi:10.1145/1391989.1391995

Timothy A. Davis (2006) Direct Methods for Sparse Linear Systems, SIAM Series “Fundamentals of Algorithms”.

Class definitions CHMfactor and dsCMatrix and function expand. Note the extra solve(*, system = . ) options in CHMfactor.

Note that chol() returns matrices (inheriting from "Matrix") whereas Cholesky() returns a "CHMfactor" object, and hence a typical user will rather use chol(A).

### Examples

data(KNex)
mtm <- with(KNex, crossprod(mm))
str(mtm@factors) # empty list()
(C1 <- Cholesky(mtm))             # uses show(<MatrixFactorization>)
str(mtm@factors) # 'sPDCholesky' (simpl)
(Cm <- Cholesky(mtm, super = TRUE))
c(C1 = isLDL(C1), Cm = isLDL(Cm))
str(mtm@factors) # 'sPDCholesky'  *and* 'SPdCholesky'
str(cm1  <- as(C1, "sparseMatrix"))
str(cmat <- as(Cm, "sparseMatrix"))# hmm: super is *less* sparse here
cm1[1:20, 1:20]

b <- matrix(c(rep(0, 711), 1), ncol = 1)
## solve(Cm, b) by default solves  Ax = b, where A = Cm'Cm (= mtm)!
## hence, the identical() check *should* work, but fails on some GOTOblas:
x <- solve(Cm, b)
stopifnot(identical(x, solve(Cm, b, system = "A")),
all.equal(x, solve(mtm, b)))

Cn <- Cholesky(mtm, perm = FALSE)# no permutation -- much worse:
sizes <- c(simple = object.size(C1),
super  = object.size(Cm),
noPerm = object.size(Cn))
## simple is 100, super= 137, noPerm= 812 :
noquote(cbind(format(100 * sizes / sizes, digits=4)))

## Visualize the sparseness:
dq <- function(ch) paste('"',ch,'"', sep="") ## dQuote(<UTF-8>) gives bad plots
image(mtm, main=paste("crossprod(mm) : Sparse", dq(class(mtm))))
image(cm1, main= paste("as(Cholesky(crossprod(mm)),\"sparseMatrix\"):",
dq(class(cm1))))

## Smaller example, with same matrix as in  help(chol) :
(mm <- Matrix(toeplitz(c(10, 0, 1, 0, 3)), sparse = TRUE)) # 5 x 5
(opts <- expand.grid(perm = c(TRUE,FALSE), LDL = c(TRUE,FALSE), super = c(FALSE,TRUE)))
rr <- lapply(seq_len(nrow(opts)), function(i)
do.call(Cholesky, c(list(A = mm), opts[i,])))
nn <- do.call(expand.grid, c(attr(opts, "out.attrs")$dimnames, stringsAsFactors=FALSE,KEEP.OUT.ATTRS=FALSE)) names(rr) <- apply(nn, 1, function(r) paste(sub("(=.).*","\\1", r), collapse=",")) str(rr, max.level=1) str(re <- lapply(rr, expand), max.level=2) ## each has a 'P' and a 'L' matrix R0 <- chol(mm, pivot=FALSE) R1 <- chol(mm, pivot=TRUE ) stopifnot(all.equal(t(R1), re[]$L),
all.equal(t(R0), re[]$L), identical(as(1:5, "pMatrix"), re[]$P), # no pivoting
TRUE)

# Version of the underlying SuiteSparse library by Tim Davis :
.SuiteSparse_version()


[Package Matrix version 1.5-3 Index]