[Rd] error message in chol() (PR#1061)

[jerome@stat.ubc.ca]

Full_Name: Jerome Asselin
Version: 1.3.0
OS: Windows 98
Submission from: (NULL) (24.77.112.193)



I am having accuracy problems involving the computation of inverse of
nonnegative definite matrices with solve(). I also have to compute the
Choleski decomposition of matrices. My numerical problems involving solve()
made me find a bug in the chol() function. Here is an example.

#Please, load the (8x8) matrix "mat" given below.
#The matrix "mat" is indeed invertible as
solve(mat)
#exists. However, the command
chol(mat)
#claims that "mat" is singular. The correct error message should be
#that "mat" is not symmetric positive definite.

#In fact, the matrix "mat" should be symmetric, but if
#we adjust for rounding errors by
chol((mat+t(mat))/2)
#we have the same error message in R, whereas S-Plus works (with
#a warning message saying that the rank is not full).

mat <- structure(c(0.587765115222886, 0.167899163888134, -1.17517413388306e-05,
  2.84014489816759e-05, -0.00314742260306333, 0.00957598738126672, 
 -0.000181283657365524, 0.000268404362197193, 0.167899163888134, 
 0.830166933867836, 2.84014489816759e-05, -6.37164529062413e-05,  
-0.00221693916036902, 0.00290212917834065, 9.0762242654496e-05, 
 -0.000148298146106774, -1.17517413388306e-05, 2.84014489816759e-05, 
 0.592681746779869, 0.155708697468761, -0.000181283657365506, 
 0.000268404362197185, -0.00203952216264033, 0.00850865103383899, 
 2.84014489816759e-05, -6.37164529062413e-05, 0.155708697468761, 
 0.856464279975445, 9.0762242654496e-05, -0.000148298146106784, 
 -0.00296804648853004, 0.00336588501549565, -0.00314742260306333,
  -0.00221693916036901, -0.000181283657365512, 9.07622426544948e-05,
  1.72579173368880, 0.238005942444999, 0.120912001671674, -0.141897803793986,
  0.00957598738126672, 0.00290212917834065, 0.000268404362197191, 
 -0.000148298146106782, 0.238005942444999, 1.56986591648628, 
-0.141897803793988,  0.213872027136185, -0.000181283657365518, 
9.07622426544981e-05,  -0.00203952216264033, -0.00296804648853004,
 0.120912001671674,  -0.141897803793987, 1.53425760353976, 1.16174845591861,
 0.000268404362197194,  -0.000148298146106774, 0.00850865103383899,
 0.00336588501549565,  -0.141897803793985, 0.213872027136185, 
1.16174845591861, 0.773168806527852 ), .Dim = c(8, 8)) 



-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._