chol | R Documentation |
Compute the Choleski factorization of a real symmetric positive-definite square matrix.
chol(x, ...) ## Default S3 method: chol(x, pivot = FALSE, LINPACK = FALSE, tol = -1, ...)
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 an error)? |
tol |
A numeric tolerance for use with |
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.
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.
This is an interface to the LAPACK routines DPOTRF
and
DPSTRF
,
LAPACK is from https://www.netlib.org/lapack/ and its guide is listed in the references.
Anderson. E. and ten others (1999)
LAPACK Users' Guide. Third Edition. SIAM.
Available on-line at
https://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.
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]) 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
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.