The Choleski Decomposition

Description

Compute the Choleski factorization of a symmetric (Hermitian), positive definite square matrix.

Usage

chol(x, pivot = FALSE)

Arguments

Details

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 may or may not occur, depending on finite precision numerical errors.

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.

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 or for non-negative definiteness. If pivot = TRUE and x is not non-negative definite then there will be no error message but a meaningless result will occur. So only use pivot = TRUE when x is non-negative definite by construction.

References

Dongarra, J. J., Bunch, J. R., Moler, C. B. and Stewart, G. W. (1978) LINPACK Users Guide. Philadelphia: SIAM Publications.

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'
all(abs(m  -  t(cm) %*% cm) < 100* .Machine$double.eps) # TRUE


x <- matrix(c(1:5, (1:5)^2), 5, 2)
m <- crossprod(x)
Q <- chol(m)
all.equal(t(Q) %*% Q, m) # TRUE

Q <- chol(m, pivot = TRUE)
pivot <- attr(Q, "pivot")
oo <- order(pivot)
all.equal(t(Q[, oo]) %*% Q[, oo], m) # TRUE
all.equal(t(Q) %*% Q, m[pivot, pivot]) # TRUE

# now for something postitive semi-definite
x <- cbind(x, x[, 1]+3*x[, 2])
m <- crossprod(x)
qr(m)$rank # is 2, as it should be

test <- try(Q <- chol(m))

(Q <- chol(m, pivot = TRUE)) # NB wrong rank here ... you have been warned!
pivot <- attr(Q, "pivot")
oo <- order(pivot)
all.equal(t(Q[, oo]) %*% Q[, oo], m) # TRUE
all.equal(t(Q) %*% Q, m[pivot, pivot]) # TRUE