The Choleski Decomposition
Compute the Choleski factorization of a real symmetric positive-definite square matrix.
chol(x, ...)"chol"(x, pivot = FALSE, LINPACK = FALSE, tol = -1, ...)
- 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.
- Should pivoting be used?
- logical. Should LINPACK be used (now ignored)?
- A numeric tolerance for use with
pivot = TRUE.
chol is generic: the description here applies to the default
Note that only the upper triangular part of
x is used, so
that $R'R = x$ when
x is symmetric.
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.
pivot = TRUE, then the Choleski decomposition of a positive
x can be computed. The rank of
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
t(Q) %*% Q equals
pivot]. See the examples.
The value of
tol is passed to LAPACK, with negative values
selecting the default tolerance of (usually)
.Machine$double.neg.eps * max(diag(x). The algorithm terminates once
the pivot is less than
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
"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
This is an interface to the LAPACK routines
DPSTRF, LAPACK is from http://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 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.
( 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