| chol {Matrix} | R Documentation |
Compute the Choleski factorization of a real symmetric positive-definite square matrix.
chol(x, ...) ## S4 method for signature 'dsCMatrix' chol(x, pivot = FALSE, ...) ## S4 method for signature 'dsparseMatrix' chol(x, pivot = FALSE, cache = TRUE, ...)
x |
a (sparse or dense) square matrix, here inheriting from class
|
pivot |
logical indicating if pivoting is to be used. Currently, this is not made use of for dense matrices. |
cache |
logical indicating if the result should be cached in
|
... |
potentially further arguments passed to methods. |
Note that these Cholesky factorizations are typically cached with
x currently, and these caches are available in
x@factors, which may be useful for the sparse case when
pivot = TRUE, where the permutation can be retrieved; see also
the examples.
However, this should not be considered part of the API and made
use of. Rather consider Cholesky() in such situations,
since chol(x, pivot=TRUE) uses the same algorithm (but not the
same return value!) as Cholesky(x, LDL=FALSE) and
chol(x) corresponds to
Cholesky(x, perm=FALSE, LDL=FALSE).
a matrix of class Cholesky,
i.e., upper triangular: R such that R'R = x (if
pivot=FALSE) or P' R'R P = x (if
pivot=TRUE and P is the corresponding permutation matrix).
Use showMethods(chol) to see all; some are worth
mentioning here:
signature(x = "dgeMatrix"): works via
"dpoMatrix", see class dpoMatrix.
signature(x = "dpoMatrix"):
Returns (and stores) the Cholesky decomposition of x, via
LAPACK routines dlacpy and dpotrf.
signature(x = "dppMatrix"):
Returns (and stores) the Cholesky decomposition via LAPACK routine
dpptrf.
signature(x = "dsCMatrix", pivot = "logical"):
Returns (and stores) the Cholesky decomposition of x. If
pivot is true, the Approximate Minimal Degree (AMD)
algorithm is used to create a reordering of the rows and columns
of x so as to reduce fill-in.
Timothy A. Davis (2006) Direct Methods for Sparse Linear Systems, SIAM Series “Fundamentals of Algorithms”.
Tim Davis (1996), An approximate minimal degree ordering algorithm, SIAM J. Matrix Analysis and Applications, 17, 4, 886–905.
The default from base, chol; for more
flexibility (but not returning a matrix!) Cholesky.
showMethods(chol, inherited = FALSE) # show different methods
sy2 <- new("dsyMatrix", Dim = as.integer(c(2,2)), x = c(14, NA,32,77))
(c2 <- chol(sy2))#-> "Cholesky" matrix
stopifnot(all.equal(c2, chol(as(sy2, "dpoMatrix")), tolerance= 1e-13))
str(c2)
## An example where chol() can't work
(sy3 <- new("dsyMatrix", Dim = as.integer(c(2,2)), x = c(14, -1, 2, -7)))
try(chol(sy3)) # error, since it is not positive definite
## A sparse example --- exemplifying 'pivot'
(mm <- toeplitz(as(c(10, 0, 1, 0, 3), "sparseVector"))) # 5 x 5
(R <- chol(mm)) ## default: pivot = FALSE
R2 <- chol(mm, pivot=FALSE)
stopifnot( identical(R, R2), all.equal(crossprod(R), mm) )
(R. <- chol(mm, pivot=TRUE))# nice band structure,
## but of course crossprod(R.) is *NOT* equal to mm
## --> see Cholesky() and its examples, for the pivot structure & factorization
stopifnot(all.equal(sqrt(det(mm)), det(R)),
all.equal(prod(diag(R)), det(R)),
all.equal(prod(diag(R.)), det(R)))
## a second, even sparser example:
(M2 <- toeplitz(as(c(1,.5, rep(0,12), -.1), "sparseVector")))
c2 <- chol(M2)
C2 <- chol(M2, pivot=TRUE)
## For the experts, check the caching of the factorizations:
ff <- M2@factors[["spdCholesky"]]
FF <- M2@factors[["sPdCholesky"]]
L1 <- as(ff, "Matrix")# pivot=FALSE: no perm.
L2 <- as(FF, "Matrix"); P2 <- as(FF, "pMatrix")
stopifnot(identical(t(L1), c2),
all.equal(t(L2), C2, tolerance=0),#-- why not identical()?
all.equal(M2, tcrossprod(L1)), # M = LL'
all.equal(M2, crossprod(crossprod(L2, P2)))# M = P'L L'P
)