In the positive semi definite case, the cross product of the pivoted choleski factor is not always (numerically) equal to the pivoted version of the original matrix, as it should be (see ?chol). Here is some example code... set.seed(0) n <- 20;p <- 6 X <- matrix(runif(n*p),n,p) X[,5] <- X[,4];X[,1] <- X[,3]+X[,4] ## make rank def R <- chol(crossprod(X),pivot=TRUE,tol=0) ## same with default tol crossprod(R)-crossprod(X[,attr(R,"pivot")]) ## should be zero ## ... largest element in result is 70.6 R[5:6,6] <- 0 crossprod(R)-crossprod(X[,attr(R,"pivot")]) ## is zero > sessionInfo() R version 3.1.0 (2014-04-10) Platform: x86_64-unknown-linux-gnu (64-bit) locale: [1] LC_CTYPE=en_GB.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_GB.UTF-8 LC_COLLATE=en_GB.UTF-8 [5] LC_MONETARY=en_GB.UTF-8 LC_MESSAGES=en_GB.UTF-8 [7] LC_PAPER=en_GB.UTF-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods base

# Here are some simple test cases. None of these should work, but Z2, Z7, and Z8 # do (Windows, 64- and 32-bit, 3.2.0) Z0 <- matrix(0,2,2); chol(Z0) Z1 <- matrix(1,2,2); chol(Z1) Z2 <- matrix(2,2,2); chol(Z2) Z3 <- matrix(3,2,2); chol(Z3) Z4 <- matrix(4,2,2); chol(Z4) Z5 <- matrix(5,2,2); chol(Z5) Z6 <- matrix(6,2,2); chol(Z6) Z7 <- matrix(7,2,2); chol(Z7) Z8 <- matrix(8,2,2); chol(Z8) Z9 <- matrix(9,2,2); chol(Z9)

Coomment #1 is not relevant to the original issue. There's nothing wrong in principle with cholesky factoring a square matrix with all elements identical to a. The result is a matrix with the top row equal to sqrt(a) and zero otherwise. The construction of chol() is that unless pivot=TRUE, you get an error if the matrix is not positive definite, but roundoff may cause that not to be detected.