chol
The Choleski Decomposition
Compute the Choleski factorization of a real symmetric positive-definite square matrix.
Usage
chol(x, …)# S3 method for default
chol(x, pivot = FALSE, LINPACK = FALSE, tol = -1, …)
Arguments
- x
an object for which a method exists. The default method applies to numeric (or logical) symmetric, positive-definite matrices.
- …
arguments to be based to or from methods.
- pivot
Should pivoting be used?
- LINPACK
logical. Should LINPACK be used (now ignored)?
- tol
A numeric tolerance for use with
pivot = TRUE.
Details
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.
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.
Unsuccessful results from the underlying LAPACK code will result in an error giving a positive error code: these can only be interpreted by detailed study of the FORTRAN code.
Value
The upper triangular factor of the Choleski decomposition, i.e., the matrix \(R\) such that \(R'R = x\) (see example).
If pivoting is used, then two additional attributes
"pivot" and "rank" are also returned.
Warning
The code does not check for symmetry.
If 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.
References
Anderson. E. and ten others (1999) LAPACK Users' Guide. Third Edition. SIAM. Available on-line at http://www.netlib.org/lapack/lug/lapack_lug.html.
Becker, R. A., Chambers, J. M. and Wilks, A. R. (1988) The New S Language. Wadsworth & Brooks/Cole.
See Also
chol2inv for its inverse (without pivoting),
backsolve for solving linear systems with upper
triangular left sides.
Examples
library(base)
# NOT RUN {
( 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
# }