Cholesky {Matrix} | R Documentation |

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.

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

`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 |

`LDL` |
logical scalar indicating if the decomposition should be
computed as LDL' where |

`super` |
logical scalar indicating if a supernodal decomposition
should be created. The alternative is a simplicial decomposition.
Default is |

`Imult` |
numeric scalar which defaults to zero. The matrix that is
decomposed is |

`...` |
further arguments passed to or from other methods. |

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' L~ D 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' R~' R~ P = (R~ P)' (R~ P),*

but *R~ P* is *not* triangular even though *R~* is.

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.

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)`

.

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), nc = 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[1], 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.attr")$dimnames, stringsAsFactors=FALSE,KEEP.OUT.ATTRS=FALSE)) names(rr) <- apply(nn, 1, function(r) paste(sub("(=.).*","\\1", r), collapse=",")) str(rr, max=1) str(re <- lapply(rr, expand), max=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[[1]]$L), all.equal(t(R0), re[[2]]$L), identical(as(1:5, "pMatrix"), re[[2]]$P), # no pivoting TRUE) # Version of the underlying SuiteSparse library by Tim Davis : .SuiteSparse_version()

[Package *Matrix* version 1.2-17 Index]