Last chance! 50% off unlimited learning
Sale ends in
chol(x, ...)
"chol"(x, pivot = FALSE, LINPACK = FALSE, tol = -1, ...)
pivot = TRUE
."pivot"
and "rank"
are also returned.
pivot = TRUE
and x
is not non-negative definite then
there will be a warning message but a meaningless result will occur.
So only use pivot = TRUE
when x
is non-negative definite
by construction.DPOTRF
and
DPSTRF
, LAPACK is from http://www.netlib.org/lapack and its guide is listed
in the references.chol
is generic: the description here applies to the default
method. Note that only the upper triangular part of x
is used, so
that $R'R = x$ when x
is symmetric.
If pivot = FALSE
and x
is not non-negative definite an
error occurs. If x
is positive semi-definite (i.e., some zero
eigenvalues) an error will also occur as a numerical tolerance is used.
If pivot = TRUE
, then the Choleski decomposition of a positive
semi-definite x
can be computed. The rank of x
is
returned as attr(Q, "rank")
, subject to numerical errors.
The pivot is returned as attr(Q, "pivot")
. It is no longer
the case that t(Q) %*% Q
equals x
. However, setting
pivot <- attr(Q, "pivot")
and oo <- order(pivot)
, it
is true that t(Q[, oo]) %*% Q[, oo]
equals x
,
or, alternatively, t(Q) %*% Q
equals x[pivot,
pivot]
. See the examples.
Pivoting with LAPACK requires LAPACK >= 3.2 and was added in R 2.15.2
(it was previously available using LINPACK). The value of tol
is passed to LAPACK, with negative values selecting the default
tolerance of (usually) nrow(x) * .Machine$double.neg.eps *
max(diag(x)
. The algorithm terminates once the pivot is less than
tol
.
Becker, R. A., Chambers, J. M. and Wilks, A. R. (1988) The New S Language. Wadsworth & Brooks/Cole.
chol2inv
for its inverse (without pivoting),
backsolve
for solving linear systems with upper
triangular left sides.( m <- matrix(c(5,1,1,3),2,2) )
( cm <- chol(m) )
t(cm) %*% cm #-- = 'm'
crossprod(cm) #-- = 'm'
# now for something positive semi-definite
x <- matrix(c(1:5, (1:5)^2), 5, 2)
x <- cbind(x, x[, 1] + 3*x[, 2])
colnames(x) <- letters[20:22]
m <- crossprod(x)
qr(m)$rank # is 2, as it should be
# chol() may fail, depending on numerical rounding:
# chol() unlike qr() does not use a tolerance.
try(chol(m))
(Q <- chol(m, pivot = TRUE))
## we can use this by
pivot <- attr(Q, "pivot")
crossprod(Q[, order(pivot)]) # recover m
## now for a non-positive-definite matrix
( m <- matrix(c(5,-5,-5,3), 2, 2) )
try(chol(m)) # fails
(Q <- chol(m, pivot = TRUE)) # warning
crossprod(Q) # not equal to m
Run the code above in your browser using DataLab