chol {base} | R Documentation |
Compute the Choleski factorization of a real symmetric positive-definite square matrix.
chol(x, ...)
## Default S3 method:
chol(x, pivot = FALSE, LINPACK = pivot, ...)
x |
an object for which a method exists. The default method applies to real symmetric, positive-definite matrices. |
... |
arguments to be based to or from methods. |
pivot |
Should pivoting be used? |
LINPACK |
logical. Should LINPACK be used in the non-pivoting case (for compatibility with R < 1.7.0)? |
chol
is generic: the description here applies to the default
method.
This is an interface to the LAPACK routine DPOTRF and the LINPACK routines DPOFA and DCHDC.
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 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.
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.
Becker, R. A., Chambers, J. M. and Wilks, A. R. (1988) The New S Language. Wadsworth & Brooks/Cole.
Dongarra, J. J., Bunch, J. R., Moler, C. B. and Stewart, G. W. (1978) LINPACK Users Guide. Philadelphia: SIAM Publications.
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.
chol2inv
for its inverse (without pivoting),
backsolve
for solving linear systems with upper
triangular left sides.
qr
, svd
for related matrix factorizations.
( 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])
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)) # NB wrong rank here - see Warning section.
## we can use this by
pivot <- attr(Q, "pivot")
oo <- order(pivot)
t(Q[, oo]) %*% Q[, oo] # recover m
## now for a non-positive-definite matrix
( m <- matrix(c(5,-5,-5,3),2,2) )
try(chol(m)) # fails
try(chol(m, LINPACK=TRUE)) # fails
(Q <- chol(m, pivot = TRUE)) # warning
crossprod(Q) # not equal to m