Description Usage Arguments Details Value References See Also Examples
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.
1 |
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 A+m*I where m is the value of |
... |
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. http://doi.acm.org/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)
.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 | 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()
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.