chol {base} | R Documentation |
The Cholesky Decomposition
Description
Compute the Cholesky factorization of a real symmetric positive-definite square matrix.
Usage
chol(x, ...)
## Default S3 method:
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 passed to or from methods. |
pivot |
logical: should pivoting be used? |
LINPACK |
logical. Defunct and gives an error. |
tol |
a numeric tolerance for use with |
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 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 Cholesky 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 Cholesky decomposition, i.e., the
matrix such that
(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.
Source
This is an interface to the LAPACK routines DPOTRF
and
DPSTRF
,
LAPACK is from https://netlib.org/lapack/ and its guide is listed in the references.
References
Anderson. E. and ten others (1999)
LAPACK Users' Guide. Third Edition. SIAM.
Available on-line at
https://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.
qr
, svd
for related matrix factorizations.
Examples
( 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