error message in chol() (PR#1061)

On Sun, 19 Aug 2001 jerome@stat.ubc.ca wrote:

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.

Yes.  Message changed.  Actually, it only uses the upper triangle,
so symmetry is assumed, as the help page says....

#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).

Well, the symmetrized matrix is not even close to positive definite, so
what does `works' mean?

eigen((mat+t(mat))/2)

$values:
[1]  2.3768888  1.8976031  1.4630238  0.9286787  0.9159388  0.5205478  0.5018393
[8] -0.1343581
 ...

L <- chol((mat+t(mat))/2)
t(L) %*% L - (mat+t(mat))/2

has one element of about 0.27.  That's just undocumented in S-PLUS.

Generally, if you want good numerical control you should be using eigen
and not chol, and if you know your matrix is symmetric, you can do better
than solve().  Have you considered the Matrix package?

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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

-- 
Brian D. Ripley,                  ripley@stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272860 (secr)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595