Bug 15853 - chol problem with positive semi definite matrix and pivoting
Summary: chol problem with positive semi definite matrix and pivoting
Status: NEW
Alias: None
Product: R
Classification: Unclassified
Component: Misc (show other bugs)
Version: R 3.1.0
Hardware: x86_64/x64/amd64 (64-bit) Linux
: P5 normal
Assignee: R-core
URL:
Depends on:
Blocks:
 
Reported: 2014-07-05 14:57 UTC by Simon Wood
Modified: 2015-06-09 21:59 UTC (History)
3 users (show)

See Also:


Attachments

Note You need to log in before you can comment on or make changes to this bug.
Description Simon Wood 2014-07-05 14:57:35 UTC
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
Comment 1 Joe V 2015-06-09 16:52:40 UTC
# 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)
Comment 2 Peter Dalgaard 2015-06-09 21:59:44 UTC
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.