I don't understand how to use the chol function in R to factor a positive semi-definite matrix. (Or I do, and there's a bug.) The documentation states:
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 ...
The following example seems to belie this description.
> x <- matrix(1, nrow=3, ncol=3)
> Q <- chol(x, pivot=TRUE)
> oo <- order(attr(Q, 'pivot'))
> t(Q[, oo]) %*% Q[, oo]
[,1] [,2] [,3]
[1,] 1 1 1
[2,] 1 1 1
[3,] 1 1 3
The result is not x. Am I using the pivot incorrectly?
For full rank input, i.e., a positive-definite matrix x, we need
Q <- chol(x, TRUE)
oo <- order(attr(Q, 'pivot'))
unpivQ <- Q[, oo]
all.equal(crossprod(unpivQ), x)
For a valid rank-deficient input, i.e., a positive semi-definite matrix x (indefinite matrix with negative eigen values are illegal, but not checked in chol), remember to zero deficient trailing diagonal block:
Q <- chol(x, TRUE)
r <- attr(Q, 'rank')
if (r < nrow(x)) Q[(r+1):nrow(x), (r+1):nrow(x)] <- 0
oo <- order(attr(Q, 'pivot'))
unpivQ <- Q[, oo]
all.equal(crossprod(unpivQ), x)
Some people call this a 'bug' of chol, but actually it is a feature of the underlying LAPACK routine dpstrf. The factorization proceeds till the first diagonal element which is below a tolerance, leaving the trailing matrix simply untouched on exit.
Thanks to Ian for the following observation:
You could use R's negative indexing in Q[-(1:r), -(1:r)] <- 0 to skip the if statement.